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ABSTRACT 


In  recent  military  conflicts,  ballistic  missiles  have  been  used  to  achieve  military 
and  psychological  objectives.  With  the  proliferation  of  weapons  of  mass  destruction 
(WMD),  the  growing  threat  of  ballistic  missiles  being  used  as  a  delivery  platform  for 
WMD  by  rogue  nations  or  militant  groups  becomes  a  concern  for  many  countries. 
Defense  against  such  threats  becomes  increasingly  important. 

There  are  different  guidance  laws  for  the  missile  interception  of  aerial  targets. 
These  include  pursuit,  proportional  navigation  (PN)  guidance  as  well  as  its  variants.  A 
new  guidance  algorithm  was  developed  by  John  A.  Lukacs  IV  and  Prof  Yakimenko  in 
2006  to  intercept  a  ballistic  missile  during  the  boost  phase  by  a  missile  interceptor.  This 
TS  guidance  algorithm  uses  the  direct  method  of  calculus  of  variations  that  maximizes 
the  kinetic  energy  transfer  from  a  surface-launched  missile  to  a  ballistic  missile  target. 

A  trade-off  study  was  conducted  by  applying  this  guidance  law  in  simulated 
ballistic  missile  interception.  This  study  examines  the  interactions  and  trade-offs  between 
the  various  critical  parameters  in  the  intercept  solution,  like  the  endgame  intercept 
geometry,  time-to-intercept  and  intercept  altitude.  It  provides  insights  into  the  feasibility 
and  limitations  of  the  TS  guidance  algorithm.  A  literature  review  of  the  drag  model  used 
in  the  algorithm  and  comparison  of  the  new  guidance  with  the  compensated  PN  guidance 
was  also  conducted.  A  new  induced  drag  model  was  developed  for  future  studies. 

The  results  verified  that  the  trajectory-shaping  guidance  is  feasible  for  the 
interception  of  ballistic  missiles  in  the  boost  phase  for  a  wide  range  of  interceptor  launch 
locations  with  respect  to  a  ballistic  missile  detection  point.  A  better  understanding  of  the 
trade-offs  between  the  key  parameters  allows  users  to  optimize  the  performance  of  this 
guidance. 
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I. 


INTRODUCTION 


A.  BACKGROUND 

In  recent  conflicts,  notably  the  first  Gulf  War  in  1992  and  the  Israeli-Lebanon 
conflict  in  2007,  the  use  of  tactical  ballistic  missiles  to  achieve  military  and  psychological 
objectives  had  alerted  political  and  military  leaders  around  the  world  of  the  immense 
threat  posed  by  such  weapons.  With  the  proliferation  of  weapons  of  mass  destruction 
(WMD),  there  has  been  a  greater  concern  that  tactical  ballistic  missiles  will  be  used  by 
rogue  nations  or  militant  groups,  as  a  delivery  platform  for  such  weapons  [1],  The 
outcome  can  be  devastating,  with  the  loss  of  many  innocent  lives  and  large-scale 
destruction.  Many  countries  have  invested  substantially  in  developing  ballistic  missile 
defense  systems  as  a  strategy  against  such  threats.  In  particular,  the  United  States 
Department  of  Defense  (DoD)  was  directed  to  develop  a  conceptual  framework  in 
2002 — the  Ballistic  Missile  Defense  System  (BMDS),  in  which  active  defense  of 
intercepting  incoming  ballistic  missiles  in  all  phases  of  its  trajectory,  is  a  major  strategy 
[2].  An  integrated  ballistic  missile  defense  system  was  subsequently  developed  based  on 
a  layered,  defense-in-depth  strategy.  In  order  to  intercept  a  ballistic  missile,  flying  at 
supersonic  speed  and  high  altitude,  the  interceptor  must  have  the  required  response  time, 
speed  and  accuracy.  The  interception  can  take  place  in  any  of  the  three  distinct  phases — 
boost,  midcourse  or  terminal,  as  illustrated  in  Figure  1 . 
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Figure  1.  Phases  of  Ballistic  Missile  Trajectory  [3] 

During  the  boost  phase,  the  ballistic  missile  rocket  motors  are  firing  in  order  to 
accelerate  the  missile  to  a  high  trajectory  altitude.  The  advantages  of  boost  phase 
interception  are  that  the  hot  and  bright  rocket  exhaust  plume  makes  detection  and 
targeting  easier,  and  decoys  cannot  be  used  during  this  phase.  The  disadvantages  lie  in 
the  geographical  siting  of  the  interceptor  system  (which  has  to  be  close  to  hostile 
territory)  and  the  short  time  to  intercept,  typically  about  180  seconds.  Figure  2  shows 
images  of  ballistic  missile  interception  by  ground-launched  intercept  missile.  In  some 
literatures,  defense  analysts  and  scientists  believed  that  interception  of  ballistic  missile  by 
surface-launched  interceptor  using  established  guidance  laws  is  not  feasible.  However, 
the  trajectory-shaping  guidance  developed  by  Lukacs  and  Prof  Yakimenko  in  2006  [4] 
shows  that  boost  phase  interception  is  possible,  with  promising  results  on  the 
effectiveness  of  such  interception,  subject  to  limitations. 
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Figure  2.  Ballistic  Missile  Interception  by  Ground-launched  Missile  [5] 

In  the  mid-course  phase,  the  ballistic  missile  is  in  space  after  the  rocket  bums  out. 
The  coast  period  through  space  before  re-entering  the  atmosphere  can  take  several 
minutes,  up  to  20  minutes  for  a  long-range  Inter-Continental  Ballistic  Missile  (ICBM). 
The  advantage  of  intercepting  a  ballistic  missile  in  this  phase  is  that  there  will  be 
sufficient  decision  and  intercept  time  as  well  as  a  greater  flexibility  in  the  geographical 
defensive  position  of  the  interceptor  system.  However,  such  interception  would  require  a 
larger  and,  hence,  heavier  interceptor  missile,  complemented  by  sophisticated  radar  and 
other  sensors  to  handle  potential  space-based  decoys. 

In  order  to  intercept  a  ballistic  missile  in  the  terminal  phase,  the  interceptor 
missile  would  have  to  do  so  after  the  ballistic  missile  re-enters  the  atmosphere.  The 
advantages  lie  in  the  requirement  for  smaller  and  lighter  interceptor  missile,  less 
sophisticated  radar  and  lower  possibility  of  decoy  as  they  are  not  likely  to  work  in  this 
phase.  The  disadvantages  are  the  very  short  reaction  time  (in  the  region  of  30  seconds  or 
less),  less  defended  geographical  coverage  and  possible  effect  of  hazardous  materials 
over  the  target  area  in  the  case  of  detonation  of  chemical,  biological  or  nuclear 
warhead(s)  mid-air. 

This  paper  focuses  on  the  interception  of  ballistic  missiles  in  the  boost  phase. 
Intercepting  a  ballistic  missile  in  its  boost  phase  is  the  ideal  solution  for  missile  defense, 
since  the  missile  is  most  vulnerable  during  this  phase  of  its  flight.  This  is  done  usually 
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over  the  launch  territory.  There  are,  however,  many  challenges  associated  with  this  phase. 
The  interceptor  will  have  to  contend  with  large  acceleration  rates,  very  short  reaction 
time  and  requires  reliable  scanning  and  tracking  [5,6].  There  are  several  systems  are 
under  development  for  conducting  boost  phase  interception,  they  include  ground-based 
missiles,  airborne  lasers  and  space-based  intercept  missiles. 

This  study  only  looks  at  surface-based  missile  interceptor.  The  missile  guidance 
law  is  one  of  the  most  important  factors  that  determines  the  feasibility  and  effectiveness 
of  the  intercept  solution.  If  intercept  is  possible,  then  the  next  step  is  the  interceptor’s 
capability  to  kill  the  target.  With  a  very  short  engagement  time,  the  interceptor  missile 
needs  a  high  speed  and  energy  to  reach  the  target.  This  imposes  a  limitation  on  the  size  of 
warhead.  Earlier  concepts  and  studies  recognized  that  that  a  simple  warhead  effect  is  not 
sufficient  to  destroy  an  ICBM,  hence  the  initiation  of  the  development  of  hit-to-kill 
technologies  [7,8].  This  concept  relies  on  relatively  smaller  high  speed  missile,  guided 
accurately  to  the  target  and  transferring  maximum  kinetic  energy  to  the  ballistic  missile 
for  a  kill.  It  suggests  some  form  of  geometry  control  during  impact. 

The  actual  impact  geometry  is  not  an  important  parameter  for  most  guidance  laws 
implemented  in  intercept  missiles  which  uses  warhead  effect  (proximity  fuze)  as  their 
main  kill  mechanism.  However,  for  a  hit-to-kill  endgame  condition  which  demands  for 
maximum  transfer  of  kinetic  energy  to  the  target  missile,  the  intercept  geometry 
commands  a  greater  attention  and  has  to  be  controlled  as  an  input  to  the  guidance  law. 

The  objective  of  this  thesis  is  to  conduct  a  trade-off  study  to  evaluate  the 
effectiveness  of  the  hit-to-kill  trajectory- shaping  guidance  through  simulation,  using  the 
guidance  algorithm  developed  by  Lukacs  and  Prof  Yakimenko  in  2006.  The  code 
developed  for  the  guidance  law  will  “generate  the  interceptor’s  entire  flight  path  in  order 
to  minimize  the  distance  traveled,  minimize  the  time  to  intercept,  and  maximize  kinetic 
energy  transfer  by  controlling  the  interception  geometry  while  providing  near-optimal 
flight  path  to  interception.  This  will  be  done  by  utilizing  the  direct  method  of  calculus  of 
variations  combined  with  inverse  dynamics  theory  to  reverse  engineer  in  real  time,  an 
optimal  flight  path  using  the  missile’s  onboard  sensors  and  computers’’.  [9] 
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B.  THESIS  ORGANIZATION 

This  thesis  consists  of  six  chapters.  Chapter  I  provides  a  brief  introduction  and 
overview  of  ballistic  missiles  interception  in  the  boost  phase  and  outlines  the  objectives 
of  this  thesis. 

Chapter  II  provides  a  literature  review  of  the  models  and  guidance  program, 
developed  by  Lukacs  and  Prof  Yakimenko  in  2006,  based  on  the  trajectory- shaping 
guidance  law.  The  trade-off  study  is  based  on  this  guidance  algorithm,  with  only  slight 
modification  in  certain  areas,  for  the  simulation.  For  completeness,  it  is  instructive  to 
provide  a  brief  overview  of  the  guidance  law,  key  characteristics  and  the  MATLAB 
program  developed  as  a  prelude  to  the  trade-off  study  and  discussion.  The  detailed 
description  of  the  theory  behind  the  guidance  algorithm  and  models  used  is  appended  in 
Appendix  E.  Much  of  the  content  in  this  chapter  is  extracted  from  the  thesis  written  by 
Lukacs  [10]. 

Chapter  III  discusses  an  alternative  guidance  law — the  well-established 
Proportional  Navigation  (PN)  guidance  that  is  used  in  many  advance  missiles  currently  in 
service.  It  provides  some  background  information  on  PN  guidance  and  makes  some 
comparison  with  the  trajectory-shaping  guidance. 

Chapter  IV  describes  the  drag  model  that  is  used  in  the  trajectory- shaping 
guidance  code  and  presents  the  need  to  incorporate  an  induced  drag  model  as  an 
enhancement  to  the  existing  model  for  better  estimation  of  drag  on  both  the  ballistic  and 
interceptor  missile. 

Chapter  V  summarizes  the  results  and  discusses  the  feasibility  of  employing  such 
guidance  in  a  real-world  scenario. 

The  Appendices  include  data  and  plots  of  the  simulation  results  obtained,  report 
on  the  author’s  participation  in  the  AIAA  Rocket  Launch  Competition,  MATLAB  code 
for  the  induced  drag  model,  detailed  description  of  the  trajectory-shaping  guidance  and 
the  complete  MATLAB  program  codes. 
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II.  SIMULATION  AND  MODELING 


A.  DESCRIPTION  OF  TRAJECTORY  SHAPING  GUIDANCE 

In  2006,  Lukacs  and  Prof  Yakimenko  [9]  developed  a  simulation  code  in 
MATLAB  as  a  demonstration  of  the  feasibility  of  intercepting  a  ballistic  missile  during 
the  boost  phase  by  a  surface-launched  interceptor  missile  using  the  Trajectory  Shaping 
guidance  law  (henceforth  termed  as  TS  guidance  in  the  report). 

The  TS  guidance  uses  the  principle  of  flight  path  optimization  from  the 
interceptor  to  the  predicted  target  position.  It  relies  on  high-order  polynomial  as  a 
reference  function  for  the  flight  path  and  uses  virtual  as  opposed  to  physical  domain  in 
optimization  process.  A  preset  thrust  history  is  used  in  the  computation  of  interceptor 
flight  path.  The  flight  path  is  derived  by  minimization  a  combined  performance  index 
including  intercept  geometry,  time-to-intercept,  penalty  on  altitude  and  dynamic 
constraints.  The  performance  index,  J  is  expressed  in  the  following  equation, 

J  =  tgo^  +  Wia(P  -  Pdesired)^  +  Pi  +  P2 

where  tgo  is  the  time-to-intercept,  p  is  the  impact  angle.  Pi  is  the  penalty  on  intercept 
altitude  and  P2  is  the  penalty  on  dynamic  constraints.  The  application  of  such  guidance  is 
particularly  suitable  for  interception  of  targets  with  relatively  fixed  predicted  trajectory, 
like  that  of  a  ballistic  missile  in  the  boost  phase. 

In  the  simulation,  the  ballistic  missile  target  was  modeled  after  the  Taepo-Dong  2 
(TPD-2)  ballistic  missile  developed  by  the  People’s  Democratic  Republic  of  Korea 
(PDRK,  also  known  as  North  Korea)  while  the  U.S.-made  SM-6  Standard  missile,  an 
upgraded  and  extended  range  version  of  the  SM-3  (shown  in  Figure  3)  was  used  as  the 
interceptor  missile.  A  3  degree-of-ffeedom  (3DoF)  mathematical  model  was  previously 
developed  and  used  to  simulate  the  trajectory  and  flight  characteristics  of  both  the 
ballistic  and  interceptor  missile  based  on  available  missile  data  from  the  open  source.  The 
intercept  path  is  continuously  calculated  onboard  the  interceptor  missile  as  a  two-point 
boundary  value  problem,  using  Direct  Methods  of  Calculus  of  Variation  to  calculate  a 
near-optimal  flight  path  and  the  control  commands  necessary  to  achieve  it  (refer  to 
Appendix  3  for  the  detailed  description  of  the  TS  guidance). 
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Figure  3.  TPD-2  ICBM  (left)  and  SM-3  Standard  Missile  (right)  [11,12] 


The  TPD-2  ballistic  missile,  developed  by  North  Korea  is  believed  to  be  an 
intercontinental  ballistic  missile  (ICBM)  capable  of  an  effective  range  of  6000  to  6500 
km,  with  a  3-stage  booster  [13].  This  put  Alaska  and  all  Asia  Pacific  countries  up  to 
northern  Australia  within  reach  of  the  TPD-2,  from  launch  sites  within  North  Korea  [see 
Figure  4],  In  order  to  have  a  feasible  solution  of  intercepting  this  missile  during  its  boost 
phase,  the  interceptor  must  be  fired  from  a  ship  or  land-based  launcher  within  range  of 
the  ballistic  missile  launch  site. 


Figure  4.  Reach  of  the  TPD-2  ICBM  Launched  from  North  Korea 

The  SM-6  (shown  in  Figure  5)  is  the  U.S.  Navy's  latest  Extended  Range  Anti-Air 
Warfare  missile  (ERAM)  that  employs  active-homing  terminal  guidance — a  key  feature 
that  allows  for  accuracy  necessary  to  intercept  an  ICBM  in  the  boost  phase.  The  SM-6 
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missile  is  the  upgraded  and  extended  range  version  of  the  SM-3,  currently  under 
development  by  Raytheon,  for  the  U.S.  Navy  as  the  next  generation  anti-ballistic  missile 
defense  system.  It  is  designed  to  be  capable  of  intercepting  ballistic  missiles  in  the  boost 
phase  if  deployed  within  range  of  the  ballistic  missile  launch  site  and  complemented  by 
suitable  sensor  systems  to  detect  the  launch.  Figure  6  shows  the  comparison  of  the  size  of 
the  TPD-2  and  SM-6. 

_ 

Figure  5.  The  SM-6  Standard  Missile  [14] 


Figure  6.  Comparison  of  the  size  of  TPD-2  and  M-6  Standard  Missile  [after  15] 

This  chapter  provides  a  brief  overview  of  the  3-dimensional  (3D)  target  model 
that  operates  in  the  Earth’s  gravitational  field,  using  available  TPD-2  data  from  the  open 
source.  A  3-D  interceptor  model,  based  on  the  SM-6  specification,  was  developed  that 
operates  in  the  Earth’s  gravitational  field,  was  also  used.  Brief  description  of  the  various 
function  files,  assumptions  made,  data  and  values  used  in  the  simulation  study  is  also 
presented  to  help  the  reader  understand  the  program  flow.  A  listing  of  the  MATLAB 
program  code  as  well  as  detailed  description  and  theory  of  the  various  functions  are 
reproduced  in  Appendices  D  and  E,  respectively,  for  completeness  and  easy  cross- 
reference. 
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B. 


SIMULATION  SOFTWARE  ARCHITECTURE 


The  general  program  flow  of  the  3DoF  simulation  eondueted  is  shown  in  Figure 
7.  The  remaining  seetions  of  this  chapter  will  briefly  describe  the  function  of  the  main 
blocks  shown  in  the  architecture,  models  used  and  the  general  flow  of  the  simulation. 
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Figure  7.  Comparison  of  the  size  of  TPD-2  and  M-6  Standard  Missile  [15] 

C.  TARGET  MODELING 

I.  Ballistic  Missile  Physical  Specifications 

Table  1  shows  the  physical  characteristics  and  specification  of  the  TPD-2  missile 


[13]. 


Overall 

Length 

32  m 

Payload 

750-1000  kg 

Range 

3500^300  km 

Stages 

2 

Stage  1 

Stage  2 

Diameter 

2.2  m 

1.335  m 

Length 

16  m 

14  m 

Launch  Weight 

-60,000  kg 

15,200  kg 

Thrust 

-103,000  kgf 

13,350  kgf 
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Thrust 

Chambers 

4,1 

Type 

LR ICBM 

Fuel  /  Oxidizer 

TM-185/AK- 

271 

TM-185/AK- 

271 

Propellant  Mass 

- 

12,912  kg 

Bum  Time 

-125  s 

110s 

Table  1.  Specifications  of  TPD-2  Ballistic  Missile 

2.  The  Ballistic  Missile  Model  Program 

The  3DoF  model  developed  by  Lukacs  and  Prof  Yakimenko  in  2006  comprises  a 
series  of  MATLAB  functions  on  an  iterative  integration  loop,  using  4  function  fdes  to 
accomplish  the  modeling.  A  brief  description  of  each  function  is  as  follows: 

1.  BRFlightS.m  -  integrates  each  time  step  to  determine  the  current  position, 
attitude,  and  aerodynamic  forces  acting  on  the  rocket/missile; 

2.  BRParamsS.m  -  determines  the  mass  of  the  ballistic  missile  and  the 
surface  reference  area; 

3.  BRDrag.m  -  determines  the  drag  coefficient; 

4.  STatmos.m  -  determines  the  properties  of  the  local  atmosphere; 

The  program  BRFlight3.m  generates  a  ballistic  flight  path  of  the  ballistic  missile 
target  based  on  the  model  developed  by  Zarchan  [16].  The  program  generates  a  3-D 
flight  path  that  is  contained  within  the  2-D  x-z  plane  shown  in  Figure  7,  where  the 
asterisks  represent  the  staging  events. 


X  10® 


Y  position  (m)  X  position  (m) 

Ballistic  Missile  Flight  Path  [Ref  ] 
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Figure  8. 


The  launch  position  and  angle  are  as  follows: 

x,=0 

To  =0 
Zq  =  Re 

0  =  85° 

where  Re  is  the  radius  of  the  Earth  (6,378,137  m)  based  on  WGS-84  system. 

The  thrust  is  assumed  fixed  for  each  stage  of  the  propulsion.  At  130  and  240 
seconds  after  launch,  the  thrust  drops  instantaneous  to  represent  the  staging  events.  At 
completion  of  the  boost  phase,  the  thrust  is  zero.  The  thrust  profde  of  the  two  stages  is 
shown  in  Figure  8. 
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Figure  9.  TPD-2  Thrust  Profile  (Boost  Phase)  [Ref] 

In  order  to  account  for  the  aerodynamic  forces,  drag  has  to  be  calculated. 
Atmospheric  conditions,  like  density  and  temperature,  are  first  determined  using  the 
function  STatmos.m.  The  missile  drag  coefficient,  CD  is  then  computed  by  the  function 
BRDragC.m.  This  is  followed  by  calculating  the  drag  force. 

The  ballistic  missile’s  mass  is  a  simple  function  of  time  (as  shown  in  Figure  9) 
and  this  is  computed  using  the  function  BRParams3.m.  The  mass  drops  sharply  during 
the  staging  events  at  130  and  240  seconds.  After  the  completion  of  the  boost  phase,  the 
mass  remains  constant  for  the  remaining  duration  of  the  flight. 
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Figure  10.  TPD-2  Ballistic  missile  Mass  (Boost  Phase  Only) 

The  final  step  of  the  ballistic  missile  model  is  to  record  all  the  data  for  the  ballistic 
missile.  This  data  will  be  called  by  the  interceptor  model  to  simulate  detection  by  the 
missile’s  onboard  sensors.  The  interceptor  missile  will  then  register  the  location  and 
velocity  of  the  target  missile  at  the  appropriate  intervals.  In  this  program,  a  60-second 
delay  is  assumed  for  the  launch  of  the  interceptor  missile  to  account  for  the  detection  and 
decision  loop. 

The  final  ballistic  fly-out  range  and  the  altitude  profile  predicted  by  the  model  is 
presented  in  Figure  10  and  11. 


Figure  1 1 .  TPD-2  Fly-out  Range 
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Time  (sec)  Time  (sec) 

Figure  12.  TPD-2  Altitude  Profile  for  Entire  Flight  (left)  and  Boost  Phase  (right) 

D.  INTERCEPTOR  MODELING 

1.  Basic  Definitions  and  Assumptions 

The  basic  specifications  of  the  SM-6  missile  are  presented  in  Table  2  [17]. 


Overall 

Length 

6.5  m 

Payload 

115  kg 

Range 

150  km 

Stages 

2 

Thrust 

Chambers 

1,1 

Type 

ERAAW 

Stage  1 

Stage  2 

Diameter 

0.53  m 

0.34  m 

Length 

1.72  m 

4.78  m 

Launch  Weight 

712  kg 

686  kg 

Thrust 

- 

- 

Fuel  /  Oxidizer 

HTPB-AP 

TP-H 1205/6 

Propellant  Mass 

468  kg 

360  kg 

Bum  Time 

6  s 

- 

Table  2.  Specifications  of  SM-6  Interceptor  Missile 

2.  Interceptor  Missile  Model  Program 

The  3DOF  model  used  in  the  program  comprises  4  MATLAB  functions  files  ran 
on  a  repeating  integration  loop.  A  brief  description  of  the  function  files  are  as  follows: 

1.  SMFlightS.m  -  integrates  each  time  step  to  determine  the  current 
position,  attitude,  and  aerodynamic  forces  acting  on  the  missile; 


14 


2.  SMParams3.m  -  determines  the  mass  of  the  interceptor  missile  and 
the  surface  reference  area; 

3.  SMDrag.m  -  determines  the  drag  coefficient; 

4.  STatmos.m  -  determines  the  properties  of  the  local  atmosphere. 

Drag  is  calculated  in  a  similar  fashion  as  that  implemented  in  the  ballistic  missile 
model.  The  SMParams.m  function  uses  the  same  methodology  as  the  BRParams.m 
function  to  calculate  the  reference  surface  area  and  mass. 

The  thrust  profile  for  the  interceptor  missile  is  also  a  two-step  curve,  just  like  the 
ballistic  missile  profile.  The  only  difference  is  the  timing  for  the  staging  event.  This  is 
shown  in  Figure  11.  The  staging  event  occurs  at  6  and  26  seconds. 

2.5 


H  1 

0.5 


Time  (sec) 

Figure  13.  SM-6  Interceptor  Thrust  Profile 

The  mass  of  the  interceptor  mission  changes  at  6  and  26  seconds  to  represent  the 
staging  events.  Figure  13  shows  the  mass  profile.  It  shows  a  distinct  drop  at  6s  and  a 
constant  mass  after  26s. 
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Figure  14.  SM-6  Interceptor  Mass  (Boost  Phase  Only) 

E.  INITIALIZATION 

For  initialization,  the  earth  geographical  data  required  was  based  on  the  WGS84 
values  presented  in  Table  3.  [18]: 


Earth's  Radius,  Re 

6,378,137  m 

Earth's  Semi-Minor  Axis,  b 

6,356,752  m 

Earth's  Flattening  (1/Ellipticity),/ 

1/298.257223563 

Earth's  Rotation  Rate  (Qze) 

7.2921 16  e-5  rad/sec 

Earth's  Gravitational  Constant,  GM 

3.986004418el4mW 

Table  3.  WGS-84  Values 

F.  COMMON  FUNCTIONS 

Two  functions  were  commonly  by  all  the  models.  They  are: 

I.  SMDragC.m  and  BRDrag.m 

The  drag  on  the  missile  is  dependent  on  two  conditions,  the  Mach  number  and  the 
phase  the  missile — ^boost  or  glide.  Mach  number  represents  the  speed  of  missile  and  the 
phase  will  determine  the  presence  or  absence  of  the  base  drag.  The  plot  of  drag 
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coefficient  with  Mach  number  of  the  missiles  is  shown  in  Figure  14  [19],  The  same 
function  is  applied  through  two  drag  files  in  the  interceptor  and  ballistic  missile  model 
respectively. 


Drag  Coefficient  by  Mach  Number  and  Boost  Phase 


Figure  15.  Drag  Coefficient  by  Mach  Number  and  Flight  Phase 

Having  obtained  the  drag  coefficient  and  computed  the  dynamic  pressure 
computed,  drag  can  be  calculated. 

2.  STatmos.m 

The  atmospheric  parameters,  which  are  dependent  on  altitude,  are  computed  by 
STatmos.m  function  file.  The  calculation  was  based  on  the  1976  standard  atmospheric 
survey,  and  includes  values  up  to  86  km  in  a  tabular  format.  Figures  15-17  show  the 
atmospheric  charts  used  by  STatmos.m  to  derive  all  the  parameters  required. 


Atmospheric  Temperature  Variation  by  Altitude 
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Figure  16. 


Atmospheric  Temperature  Variation  by  Altitude 


Atmospheric  Density  Variation  by  Altitude 


Figure  17.  Atmospheric  Density  by  Altitude 


Atmospheric  Pressure  Variation  by  Altitude 


Figure  18.  Atmospheric  Pressure  Variation  by  Altitude 

G.  THE  FLIGHT  PROGRAM 

Once  the  previous  iteration  (or  the  initialization)  has  been  integrated,  it  is  returned 
to  the  program  as  the  current  values.  The  program  then  uses  these  values  to  calculate  all 
the  descriptive  values  of  the  system,  apply  the  corrective  time-  and  position-  dependant 
factors,  and  calculate  the  derivatives  for  the  next  iteration. 

The  BRParams.m  and  SMParams.m  functions  were  called  to  determine  the 
reference  areas  for  the  control  surfaces  and  missile  plan  form  areas.  The  program  then 
calls  the  function  ACoeff.m  to  determine  several  necessary  aerodynamic  coefficients. 
The  function  inputs  are  angle  of  attack,  altitude,  and  pitch  control  surface  deflection. 
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The  respective  missile  model  will  call  the  function  SMDrag.m  and  BRDrag.m  to 
determine  the  value  of  the  drag  coefficient.  Drag  and  thrust  are  then  calculated  and  the 
execution  is  returned  to  the  main  program. 

The  trajectory  of  the  interceptor  missile  is  generated  by  the  function 
SMTrajectory.m  by  applying  the  TS  guidance  law  through  the  SMGuidance.m  function. 
This  is  an  iterative  process,  starting  with  a  ‘guess’  of  the  interceptor  final  states  and 
subsequently  performing  trajectory  optimization  by  calling  sub-routine 
SMGuidanceCost.m  to  minimize  the  cost  and  penalty  of  each  iterative  flight  path 
generated.  The  optimized  flight  path  is  then  returned  to  SMTrajectory  for  execution 
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III.  COMPARISON  OF  COMPENSATED  PROPORTIONAL 
NAVIGATION  WITH  TRAJECTORY-SHAPING  GUIDANCE 


In  order  to  understand  the  guidance  laws,  it  is  useful  to  understand  what 
“guidance”  is.  Guidance  is  the  logic  that  issues  steering  commands  to  the  vehicle  to 
accomplish  certain  flight  objectives’  [20],  In  the  case  of  a  missile  and  its  target,  it  is  the 
algorithm  the  missile  uses  to  guide  itself  in  an  intercept  path  with  the  target,  through 
autopilot  commands  given  to  its  control  surfaces,  like  fins  or  wings.  In  current  interceptor 
missiles  being  fielded  or  deployed,  a  number  of  guidance  laws  were  used.  The  type  of 
guidance  law  implemented  depends  on  the  mission  and  the  type  of  targets  the  missile  is 
designed  for.  Examples  of  better-known  guidance  laws  are  the  Beam  Rider,  Pure  Pursuit, 
Proportional  Navigation  and  its  variants. 

A.  PROPORTIONAL  NAVIGATION  GUIDANCE 

PN  guidance  or  its  variants,  like  augmented  or  compensated  PN,  are  the  most 
common  guidance  law  implemented  in  modem  missiles.  PN  guidance  has  proven  itself  in 
many  operations  to  be  effective  against  maneuvering  target,  especially  in  an  air-to-air 
duel  or  surface-to-air  interception.  It  can  be  found  in  many  air-to-air,  air-to-surface  and 
surface-to-air  applications.  In  proportional  navigation,  the  fundamental  principle  lies  in 
the  rate  of  change  of  the  missile  heading  being  kept  proportional  to  the  rate  of  rotation  of 
the  line-of-sight  (LOS)  from  the  missile  to  the  target.  The  guidance  system  points  the 
differential  velocity  vector  at  the  target  [7].  Figure  18  illustrates  the  geometry  of 
interception  using  the  PN  guidance. 
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Figure  19.  Intercept  Geometry  for  PN  Guidance  [after  21] 

In  order  for  interception  to  occur,  the  missile  heading  angle,  'P  and  the  LOS  angle, 
X,  must  be  kept  constant.  If  the  target  increases  speed,  then  X  will  increase  and  the 
interceptor  must  correspondingly  increase  to  maintain  its  interception  course.  The  rate 
of  change  of  both  angles  must  be  proportional,  given  by  the  expression  [16], 

^  =  nA  (III.A.l) 

where  N  is  the  proportionality  constant,  which  depends  on  the  target  lead  factor  of  the 
interceptor  (usually  between  the  values  of  3-5  depending  on  the  interceptor  missile’s 
maneuverability).  The  interceptor  will  maneuver  until  1  =  0  (no  more  changes  in  LOS 
rate  with  the  target),  at  which  point  'P  =  0  .  The  interceptor  continuously  adjusts  itself  to 
maintain  the  above  condition  till  interception  (or  a  miss)  occurs. 

During  the  intercept  process,  the  missile  seeker  continuously  points  itself  at  the 
target  to  derive  the  LOS  rate  and  provides  a  signal  to  the  guidance  computer.  The 
guidance  computer  in  turn  will  convert  the  feedback  into  an  acceleration  command  to  the 
control  surfaces  to  physically  steer  the  missile.  The  normal  acceleration  to  the  missile’s 
velocity  vector  is  given  by 

a  =  (IILA.2) 

which  can  be  correlated  to  the  LOS  rate 

a  =  NV^A  (IILA.3) 
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This  is  applicable  to  a  2-D  case.  Zipfel  [20]  develops  a  3-D  dimensional  PN  guidance  law 
with  the  expression 

a  =  NV^QP^u^-g  (III.A.4) 

where  is  the  cross  product  of  the  LOS  frame  with  respect  to  the  inertia  Earth  frame, 
is  the  unit  vector  of  ,  and  g  is  the  added  gravity  bias,  which  counteracts  the  sagging 

tendency  of  the  trajectory  under  seeker  control.  This  form  of  guidance  is  termed  by  Zipfel 
as  Pure  PN. 

B.  COMPENSATED  PROPORTIONAL  NAVIGATION 

Pure  PN,  however,  is  not  commonly  implemented  in  its  original  form.  The  thrust 
generated  by  the  interceptor’s  propulsion  creates  a  parasitic  acceleration  in  the  LOS  angle 
that  has  to  be  compensated  in  the  autopilot  command.  Guidance  in  this  form  is  termed 
compensated  PN.  The  parasitic  acceleration  projected  onto  the  LOS  plane  can  be 
expressed  as  follows: 

(iii-B.i) 

It  is  then  subtracted  from  the  PN  command  in  equation  (3.B.4)  to  obtain  the  compensated 
command  [20] 

a  =  NV°  (III.B.2) 

where  the  rotation  matrix  of  the  LOS  coordinates  with  respect  to  the  wind  frame  is 
defined  by  the  azimuth  and  elevation  angles  from  the  LOS  vector. 

C.  DISADVANTAGES  OF  COMPENSATED  PN  GUIDANCE  FOR  BOOST 
PHASE  INTERCEPTION 

Three  disadvantages  of  the  compensated  PN  for  the  boost  phase  interception  of 
ballistic  missiles  were  identified  by  Lukacs  and  Prof  Yakimenko.  The  first  disadvantage 
is  the  inherent  control  system  time  constant  of  the  PN  guidance,  which  utilizes  current 
target  information  in  a  homing  guidance  loop  with  feedback  control.  Such  control  loops 
will  inevitably  result  in  a  finite  time  constant  which  can  result  in  considerable  miss 
distance.  The  second  is  the  disregard  of  the  end-game  environment.  The  PN  guidance  law 
is  focused  on  maintaining  the  LOS  rate  constant  continuously  based  on  current  missile 
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and  target  parameters.  However,  the  guidance  law  does  not  attempt  to  optimize  missile 
speed  and  altitude,  which  is  critical  for  the  endgame  condition,  to  ensure  that  the  missile 
interceptor  has  sufficient  acceleration  left  to  intercept  the  target.  Rather,  PN  guidance 
adopts  the  most  direct  and  minimum  acceleration  collision  path  with  the  hope  that  there  is 
little  target  maneuver  and  hence  sufficient  acceleration  left  in  the  interceptor  for  the 
endgame.  Lastly,  the  intercept  geometry  is  not  controlled  for  the  case  of  PN  guidance, 
which  is  left  to  the  relative  speeds  of  the  interceptor  and  target  as  well  as  the  endgame 
maneuvers  involved  Controlling  the  intercept  geometry  allows  for  the  maximization  of 
kinetic  energy  transfer  to  the  target,  which  is  critical  for  boost  phase  intercept  of  ballistic 
missile  since  a  lightweight  (very  small  or  no  warhead),  high  speed  missile  with  sufficient 
range  (most  of  the  weight  goes  to  the  fuel)  is  required. 

D.  ADVANTAGES  OF  TRAJECTORY  SHAPING  GUIDANCE 

In  TS  guidance,  the  derivation  of  the  intercept  solution  -  the  optimized  interceptor 
flight  path,  results  in  a  set  of  functions,  which  occurs  through  a  Cost  Function  and  a 
Penalty  function,  which  must  be  minimized  through  multiple  iterations.  The  three 
parameters  represented  in  the  Cost  Function  are  the  length  of  the  virtual  arc  (proportional 
to  the  flight  path  distance),  the  time  to  intercept  and  the  impact  angle  of  the  final  intercept 
(angular  deviation  from  desired  90°  impact).  Each  of  the  parameters  needs  to  be  carefully 
weighted  to  correctly  reflect  the  desired  condition  of  the  intercept,  and  affects  the  Cost 
Function  value  as  a  whole. 

The  penalty  function  consists  of  the  maximum  acceleration  in  the  y-  and  z- 
directions.  They  represent  the  ‘penalty’  to  pay  when  certain  physical  limitations  are 
reached  or  violated.  The  penalty  function  has  been  set  to  the  certain  values,  which  are 
dependent  on  speed  and  altitude  of  the  interceptor,  to  represent  the  physical  limits  beyond 
which  the  intercept  solution  will  not  be  feasible. 

The  flight  path  optimization  is  done  by  solving  a  non-linear  programming 
problem  numerically  real-time  and  once  the  minimum  function  is  obtained,  the  algorithm 
will  return  the  required  control  time  history  to  the  missile  guidance  system,  which  can 
then  execute  the  commands  and  fly  the  derived  flight  path.  Since  the  missile  system  can 
be  programmed  with  sufficient  data  to  compensate  for  its  control  system  time  constant, 
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the  system  lag  can  be  effectively  negated,  thus  eliminating  a  source  of  error.  The 
guidance  system  can  be  updated  every  few  seconds  to  increase  the  accuracy  of  the 
intercept. 

The  main  advantages  of  this  guidance  are  three-fold — (1)  the  relatively  short 
computation  time  to  iterate  and  converge  to  an  optimized  solution,  (2)  the  cost  and 
penalty  functions  are  scalable  as  desired  to  fit  the  mission  profile  and  (3)  elimination  of 
the  control  system  time  constant.  For  a  boost  phase  intercept  mission,  the  TS  guidance  is 
able  to  address  the  disadvantages  of  the  computed  PN  guidance  and  offers  greater 
flexibility  in  being  able  to  ‘customize’  the  guidance  to  improve  its  performance  to  meet 
the  different  operational  demands. 
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IV.  INDUCED  DRAG  MODEL 


A.  AERODYNAMIC  DRAG 

The  modeling  of  aerodynamic  forces  acting  on  a  missile  is  a  necessary  step  in  the 
simulation  of  the  missile  interception  of  an  aerial  target.  Drag,  which  represents  the  total 
sum  of  forces  that  opposes  the  thrust  provided  by  the  missile  propulsion  system,  is  a 
‘major  design  parameter  in  satisfying  the  flight  range  requirement  of  tactical  missiles’ 
[22],  This  is  especially  so  for  supersonic  missile.  It  is  part  of  the  guidance  algorithm  to 
incorporate  a  drag  model  to  account  for  the  effect  of  drag  on  the  missile,  which  affect 
flight  performance  in  all  phases.  A  suitable  drag  model  will  enhance  the  accuracy  of  the 
simulation  result  and  help  to  provide  a  more  realistic  prediction  of  the  flight  characteristic 
and  performance. 

B.  COMPONENTS  OF  DRAG 

In  the  earth  atmosphere,  any  flying  object  will  experience  drag.  It  is  the  force  that 
opposes  thrust  (or  the  forward  motion)  of  aircraft,  missiles  or  other  flying  platforms.  The 
overall  drag  on  a  body  is  made  up  of  several  components — friction,  base,  wave  (or 
pressure)  and  induced  drag.  Friction  drag  refers  to  the  drag  generated  by  the  skin  friction 
between  tangential  air  flow  and  the  moving  body.  It  is  primarily  a  function  of  the  body 
fineness  ratio  or  length-to-diameter  ratio  of  the  missile  body.  Base  drag  accounts  for  the 
pressure  difference  of  air  flow  between  the  nose  and  the  base  of  the  body.  It  is 
particularly  significant  during  coast  or  glide  phase  of  a  missile  since  the  propulsion 
system  is  no  longer  generating  thrust.  Wave  drag  is  caused  by  the  effect  of  normal 
pressure  of  air  and  largely  dependent  on  the  shape  of  the  nose.  The  nose  fineness  ratio 
gives  a  measure  of  the  nose  shape,  which  affects  the  wave  drag.  A  sharp  nose  will  create 
less  drag  than  a  blunt  one.  Last  but  not  least  is  induced  drag.  This  is  the  component  of 
drag  that  is  ‘induced’  by  the  lift  generated  for  level  flight  or  climb.  It  is  directly 
proportional  to  the  lift.  Hence  induced  drag  is  present  in  all  flying  platform  as  long  as  lift 
is  generated  by  the  control  surfaces  and  in  some  condition,  is  a  key  contributor  to  the 
overall  drag. 
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In  mathematical  terms,  drag  is  a  function  of  the  drag  coefficient,  dynamic 
pressure  and  reference  area,  given  by  this  equation 


(IV.B.l) 


The  drag  coefficient,  Cd  is  usually  derived  from  empirical  equation  as  a  function 
of  the  speed  of  sound  or  Mach  Number  (M)  in  the  operating  regime,  the  angle  of  attack, 
a,  which  is  the  angle  of  inclination  between  the  missile  body  axis  and  its  velocity  vector. 
In  the  case  of  a  missile,  the  reference  area  is  commonly  taken  as  the  cross-sectional  area. 

The  atmospheric  dynamic  pressure,  q  is  given  by 


(IV.B.l) 


where  p  is  the  density  of  air  and  V  is  the  speed  of  the  missile. 

C.  INDUCED  DRAG  POLAR 

As  briefly  described  in  Chapter  2,  for  the  new  TS  guidance,  a  standard  drag  model 
was  used  whereby  the  equation  (IV.B.l)  was  used  to  calculate  the  drag.  Noticed  it  was 
mentioned  that  Cd  is  usually  obtained  empirically  and  different  missile  design  is  likely  to 
have  different  drag  characteristics.  The  accuracy  of  the  drag  prediction  can  be  enhanced 
with  the  use  of  an  induced  drag  model.  This  model  predicts  drag  by  involving  the  lift 
coefficient.  Cl  which  describes  the  lift  generated  by  the  missile.  The  lift  force,  L  is 
normal  to  the  missile  velocity  vector  and  hence  also  normal  to  drag,  D  (which  is  opposite 
to  the  velocity  vector).  The  lift  and  drag  coefficient  can  be  expressed  as 


L 


(IV.C.l) 


And 


P 


(IV.C.2) 


Both  coefficients  can  be  assumed  to  be  functions  of  the  following  parameters: 
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Cl,  Cd  =f  (Mach  number,  angel  of  attack,  power  on/off,  shape)  [20] 

The  Mach  number  can  have  a  great  effect  on  the  coefficients,  especially  in  the 
transonic  and  supersonic  regime.  Of  note,  at  the  transonic  regime,  there  is  a  sudden  steep 
rise  in  drag  before  the  speed  crosses  the  sonic  line,  a  phenomenon  commonly  known  as 
transonic  drag  rise  (shown  in  Figure  19).  The  main  effect,  however,  is  caused  by  the 
angle  of  attack.  A  slight  change  can  result  in  significant  change  in  lift  coefficient. 


When  the  lift  and  drag  coefficient  are  plotted  against  each  other,  a  near  parabolic 
curve  emerges,  for  any  given  Mach  number.  The  relationship  can  be  expressed  as: 

Cd  =  Cdo  +  k(CL  -Ciof  (IV.C.3) 

Graphically,  the  plot  is  shown  in  Figure  19  and  is  known  as  the  parabolic  drag 

polar: 


Q 


Parabolic  Drag  Polar  for  Asymmetric  Body 
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Figure  21. 


In  the  drag  polar  expression,  Cdo  is  the  zero-lift  drag  coefficient,  a  term  that 
accounts  for  the  other  drag  components  when  lift  is  zero  while  Clo  represents  the  lift 
coefficient  at  the  lowest  drag.  The  coefficient  k  is  referred  to  as  the  induced  drag 
coefficient. 

If  minimum  drag  occurs  at  zero  lift,  then  Clo  is  z:ero  and  the  parabola  is 
symmetrical  about  the  Cd  axis  and  centered  at  Clo  =  0.  This  case  is  representative  of  a 
missile,  which  is  axis-symmetric  in  nature.  The  different  points  on  the  drag  polar 
represent  the  different  angles  of  attack,  which  give  rise  to  different  lift  and  drag 
coefficient.  The  parabolic  drag  polar  of  a  missile  is  as  shown  in  Figure  20: 


Figure  22.  Parabolic  Drag  Polar  for  Symmetric  Body 

In  order  to  implement  the  induced  drag  model  in  the  algorithm  to  calculate  drag 
(D),  the  set  of  Cl-Cd  data  for  the  interceptor  and  target  at  various  Mach  number  was 
obtained  from  open  sources  and  set  up  in  the  form  of  a  look-up  table.  The  data  is  then 
curve  fitted  using  a  polynomial  (parabolic)  function.  A  plot  of  the  family  of  curves  for 
different  Mach  number  using  Excel  plot  is  shown  in  Figure  21 : 


30 


Figure  23.  Parabolic  Drag  Polar  for  Different  Mach  Number  (Symmetric  Body) 


The  values  of  Cdo  and  k  can  be  derived  from  the  coefficient  of  the  polynomial 
function.  These  values  can  be  substituted  into  equation  (IV.C.3).  Since  Clo  is  zero  for  an 
axis-symmetric  missile,  the  only  outstanding  variable  that  has  to  be  computed  is  the  lift 
coefficient,  Cl.  For  a  particular  iteration,  Cl  can  be  expressed  as  a  function  of  load  factor, 
Hz.  The  load  factor  of  a  missile  is  defined  as: 

71^  -  ^  (1V.C.4) 


Where  L  is  the  lift  and  W  is  the  weight  of  the  missile.  From  equation  (IV.C.l),  we 
can  express  L  in  terms  of  Cl- 

L  =  CiqSref  (1V.C.5) 

Therefore,  substituting  equations  (1V.C.5)  into  (1V.C.4),  Cl  can  be  computed: 


■  rtf 


(1V.C.6) 


nr 


The  overall  drag  coefficient  can  then  be  calculated  using  equation  (1V.C.3),  with 
values  of  Cdo.  k  and  Cl. 

A  MATLAB  program  was  developed  to  implement  the  induced  drag  model  in  the 
TS  guidance  algorithm.  The  code  is  appended  in  Appendix  C.  Minor  modifications  have 
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to  be  made  in  the  other  sub-routines  in  order  to  integrate  the  new  drag  model  into  the 
algorithm.  It  is  intended  to  replace  the  generic  drag  model  currently.  Though  the  new 
model  has  been  integrated  and  initially  tested  to  be  working,  results  has  not  been 
consistent.  Further  testing  using  the  new  code  and  more  simulation  scenario  to  fine-tune 
the  program  can  be  conducted  in  future  studies. 
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V.  RESULTS  AND  DISCUSSIONS 


A  total  of  270  simulations  were  conducted  in  this  trade-off  study.  The  main 
objectives  of  the  simulation  study  are:  (1)  to  verify  that  the  TS  guidance  can  be  used  in  a 
hit-to-kill  interceptor  missile  to  effectively  engage  a  ballistic  missile  in  the  boost  phase, 
(2)  to  predict  the  region  of  the  interceptor  launch  position  in  relations  to  the  ballistic 
missile  launch  point  which  can  result  in  a  feasible  intercept  solution  and  (3)  the  impact  on 
the  time-to-intercept,  altitude  and  impact  angle  deviation  by  varying  the  intercept 
geometry  coefficient,  Wia  of  the  cost  function.  In  this  study,  simulations  were  run  for 
various  interceptor  missile  (simulated  by  the  surface-launch  SM-6  missile)  launch 
positions,  for  different  northing  and  easting  against  a  ballistic  missile  (simulated  by  the 
TPD-2  ICBM)  launched  from  a  fixed  site  (arbitrarily  fixed  at  point  [0,0,Re]  in  the  earth 
coordinate  system,  whereby  Re  is  the  radius  of  the  earth).  The  ballistic  missile  is 
launched  towards  a  target  in  the  north  direction. 

A.  FEASIBILITY  OF  THE  TRAJECTORY  SHAPING  GUIDANCE 

The  simulation  results  showed  that  the  trajectory  guidance  algorithm  worked  well 
for  the  intercept  scenario  ran  and  are  within  expectation.  The  guidance  provides  feasible 
solutions  for  the  interception  of  the  TPD-2  ballistic  missile  by  the  surface  launched  SM-6 
missile.  Out  of  the  270  simulation  runs,  successful  intercept  was  achieved  for  180  runs. 
The  remaining  ‘non-feasible’  intercepts  were  due  to  the  inability  of  the  guidance 
algorithm  to  converge  to  an  optimized  problem  within  the  maximum  number  of  iterations 
allowed.  The  maximum  iterations  condition  was  set  to  act  as  the  upper  bound  and  ‘time¬ 
out’  condition,  reflection  the  real  world  situation  whereby  the  intercept  scenario  does  not 
offer  a  solution  fast  enough  to  ensure  a  successful  intercept.  This  is  realistic  as  the 
interceptor  missile  has  to  compute  the  optimized  flight  path  after  launch  and  be  guided  to 
the  predicted  intercept  point  within  a  very  small  time  window.  Such  scenario  occurs 
when  the  interceptor  system  is  ‘out  of  range’  or  within  the  ‘blind  range’  to  effectively 
intercept  the  ballistic  missile  within  the  timeframe  of  60  to  180  seconds  after  launch. 
Table  22  to  24  presents  a  summary  of  simulation  results,  showing  the  number  of 
iterations  the  algorithm  takes  to  find  the  optimized  intercept  flight  path  for  values  of  Wia 


33 


set  to  0,  100  and  1000  respectively.  The  region  in  green  represents  feasible  intercepts 
predicted  by  the  TS  algorithm  and  the  numbers  refers  to  the  number  of  iterations  required 
to  find  the  optimized  intercept  flight  path.  The  letter  ‘EX’  in  the  red  region  represented 
the  non-feasible  region  whereby  there  is  no  feasible  intercept  solution  if  the  interceptor 
missile  is  launched  from  those  positions. 


Table  4.  No  of  Iterations  Required  for  Feasible  Interceptor  Solution  (WIA  =  0) 


Table  5.  No  of  Iterations  Required  for  Feasible  Interceptor  Solution  (WIA  =  100) 
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Table  6.  No  of  Iterations  Required  for  Feasible  Interceptor  Solution  (Wia  =  1000) 


The  results  obtained  successfully  verified  that  the  TS  guidance  algorithm  does 
provide  feasible  solution  to  the  boost  phase  intercept  problem.  The  new  guidance  is 
suitable  for  scenario  involving  a  high  speed  hit-to-kill  missile  interceptor  and  a  fixed 
trajectory  ballistic  missile  target.  By  comparing  the  results  for  the  different  intercept 
geometry  coefficient  values,  it  is  observed  that  the  number  of  iterations  is  increased 
across  the  board  when  Wia  is  increased  from  0  to  1000.  This  means  that  when  a  higher 
weightage  is  placed  on  intercept  geometry  in  the  intercept  solution,  it  increases  the 
computing  resources  required  and  the  optimization  process  takes  a  longer  time  to 
converge. 


B.  REGION  OF  POSSIBLE  DEPLOYMENT  POSITION  OF  THE 
INTERCEPTOR  MISSILE 

The  simulation  study  was  able  to  provide  a  mapping  of  the  feasible  intercept 
region,  in  terms  of  the  interceptor  missile  launch  position  with  respect  to  that  of  the 
ballistic  missile.  As  the  scenario  is  symmetrical  about  the  north-south  direction,  only  half 
of  the  360-degree  coverage  needs  to  be  tested.  The  other  half  is  identical  and  a  mirror 
image  of  the  results  obtained.  The  simulation  results  obtained  and  plot  is  appended  in 
Appendix  1.  Figure  22  shows  the  feasible  intercept  region  in  green  and  non-feasible 
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region  in  red,  for  the  scenario  tested.  The  mapping  was  done  by  running  simulations  and 
collecting  intercept  data  at  intervals  of  20  km  in  the  northing  and  easting  direction. 


-140-120-100  -80  -60  -40  -20  0  20  40  60  80  100  120  140 


Easting  (km) 


Figure  24.  Region  of  Interceptor  Position  for  Interception  of  Ballistic  Missile 

It  should  be  noted  that  the  shape  of  the  boundary  has  no  special  significance 
except  for  the  fact  it  is  purely  a  2-D  plotting  limitation  of  the  software  used.  A  smooth 
curve  should  be  the  correct  representation  for  the  boundary  of  the  respective  regions. 
From  the  mapping,  several  interesting  observations  were  noted.  As  expected,  the  plot 
shows  an  off-set  in  the  feasible  intercept  region  with  respect  to  the  ballistic  missile  launch 
point.  This  represents  a  distinct  difference  in  the  intercept  range  between  front-  or  rear- 
quarter  engagement  scenarios.  One  would  expect  the  rear-sector  engagement,  represented 
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by  the  intercept  launch  points  south  of  ballistic  missile  launch  point  at  [0,0],  to  have  a 
feasible  region  with  shorter  ranges.  This  is  consistent  with  a  typical  ‘tail-chase’  intercept 
scenario,  which  requires  the  interceptor  to  be  closer  to  the  target.  The  front  sector 
engagement  is  characterized  by  longer  ranges  as  depicted  by  the  larger  area  in  the  plot. 

There  is  a  central  non-feasible  region  (in  red)  and  for  near  ranges  close  to  the 
ballistic  missile  launch  point.  This  indicated  that  certain  aerodynamic  limits  were 
exceeded,  like  the  g-  and  lateral  acceleration  limits,  which  resulted  in  a  high  penalty 
function  value  and  divergence  during  the  optimization  process.  The  latter  represents  some 
kind  of  ‘blind-range’  whereby  intercept  is  not  feasible  due  to  the  target  being  too  close 
for  the  interceptor  to  effectively  maneuver  and  to  achieve  an  intercept. 

The  usefulness  of  mapping  the  region  of  feasible  intercept  is  that  it  provides 
critical  information  for  the  operational  planners  on  the  possible  region  to  deploy  an 
interceptor  system  (using  the  TS  guidance)  given  a  ballistic  missile  threat.  The  creation  of 
such  plots  can  be  automated  and  done  in  a  short  time  using  computer  program.  The  input 
requirements  are  the  specifications  of  the  interceptor  missile  and  intelligence  on  the 
specifications/characteristics  of  the  target  missile.  The  accuracy  of  such  prediction 
software  would  depend  on  the  number  of  simulation  run  (proportional  to  the  ‘resolution’ 
or  distance  interval  required)  and  accuracy  of  the  input  data  and  model  used.  There  is 
scope  to  develop  a  computer  program  to  generate  intercept  region  plots  for  operational 
planning  purposes  in  future  studies. 

C.  EFFECT  OF  VARYING  THE  INTERCEPT  GEOMETRY  COEFFICIENT 

One  of  the  main  objectives  of  the  trajectory-shaping  guidance  is  to  effect 
maximum  transfer  of  kinetic  energy  from  the  interceptor  missile  to  the  target  missile  in  a 
hit-to-kill  interception  scenario.  This  can  be  achieved  when  the  impact  angle  onto  the 
target  missile  is  at  a  right  angle  to  the  ballistic  missile  body.  The  TS  guidance  algorithm 
allows  the  missile  designer  to  adjust  the  weightage  placed  on  the  intercept  geometry, 
amongst  other  critical  parameters,  by  varying  the  value  of  Wia,  which  is  to  be  optimized 
(minimized).  The  trade-off  study  examines  the  effect  on  time-to-go,  intercept  altitude  and 
impact  angle  deviation  due  to  changes  in  Wia. 
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In  this  study,  simulations  were  run  for  three  values  of  Wia,  namely  0,  100  and 
1000.  By  having  Wia=  0  implies  that  the  optimized  flight  path  does  not  take  into  account 
the  end-game  intercept  geometry.  In  the  case  of  Wia  =  1000,  a  higher  weightage  is  placed 
on  impact  angle  at  intercept  and  the  flight  path  will  be  optimized  to  achieve  an  impact 
angle  as  close  to  90  degrees  as  possible.  When  a  higher  ‘premium’  is  placed  on  the 
intercept  geometry,  the  algorithm  will  attempt  it  at  the  expense  of  other  performance 
parameter.  Here,  there  is  always  some  form  of  trade-off  that  a  missile  designer  has  to  be 
conscious  of 

A  comparison  of  the  time-to-go  or  time-to-intercept,  intercept  altitude  and  impact 
angle  deviation  for  different  values  of  Wia  provides  valuable  insights  into  the  effects  of 
Wia  on  the  missile  performance  under  the  same  intercept  scenario.  Table  7  presents  a 
summary  of  the  comparison  results. 


Time-to-Intercept  (s) 

Intercept  Altitude  (km) 

Impact  Angle  Dev  (°) 

WIA 

Min 

Max 

Ave 

Min 

Max 

Ave 

Min 

Max 

Ave 

0 

17.6 

61.0 

38.67 

22.3 

60.0 

38.77 

2.4 

30.3 

18.45 

100 

17.6 

61.0 

38.48 

22.3 

60.0 

38.77 

2.4 

30.3 

18.45 

1000 

17.8 

60.2 

39.00 

22.4 

59.1 

38.94 

0.0 

30.0 

10.02 

Table  7.  Time-To-Intercept,  Altitude  And  Impact  Angle  Deviation 
For  Different  Intercept  Geometry  Coefficient,  Wia 

Generally,  it  can  be  seen  that  there  is  no  significant  different  between  the  results 
obtained  for  Wia  being  set  at  0  and  100,  for  the  time-to-intercept,  intercept  altitude  and 
impact  angle  deviation.  However,  some  different  is  observed  when  the  weightage 
coefficient  is  increased  by  a  factor  of  1000.  The  most  significant  difference  is  observed  in 
the  impact  angle  deviation.  With  a  high  weightage  placed  on  intercept  geometry  (Wia  = 
1000),  the  minimum  impact  angle  deviation  is  0,  and  the  average  over  all  the  feasible 
intercepts  reduces  from  about  18°  to  10°'  This  is  within  expectation  and  verified  that  the 
algorithm  indeed  attempted  to  minimize  the  impact  angle  deviation  when  the  Wia  is 
increased.  Another  important  observation  is  that  by  enhancing  the  intercept  geometry 
weightage,  it  does  not  affect  time-to-intercept  and  intercept  altitude.  By  correlating  the 
observations  made  in  the  previous  section,  it  can  be  seen  that  when  Wia  is  increased,  the 
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number  of  iterations  needed  for  convergence  also  goes  up.  However,  the  overall  time 
taken  for  the  intercept  generally  does  not  increase  significantly.  This  seems  possible  if 
there  is  some  trade-off  or  interaction  with  the  other  parameter(s)  which  is  not  investigated 
in  this  study.  It  can  also  be  seen  that  the  area  of  intercept  decreases  when  Wia  is  changed 
from  0  to  1000.  A  possible  reason  is  that  the  larger  number  of  iterations  required  for 
higher  intercept  geometry  coefficient  may  have  resulted  in  the  further  range  intercepts 
exceeding  the  maximum  iteration  limit  imposed  by  the  guidance.  This  will  cause  the 
longer  range  region  to  be  marked  as  non-feasible.  However,  this  clearly  shows  the 
‘penalty’  for  imposing  a  placing  a  higher  weightage  on  intercept  geometry.  Figure  23 
provides  an  illustration  of  the  typical  intercept  geometry  at  Wia  =  1000.  It  can  be  seen 
that  the  flight  path  is  optimized  for  a  near  90°  impact  to  the  ballistic  missile  trajectory. 


Figure  25.  Typical  Intercept  Geometry 

It  can  be  expected  that  when  the  weightage  for  intercept  geometry  is  increased, 
there  is  greater  constraints  imposed  on  the  interceptor  system  to  optimized  impact  angle 
amongst  other  factors.  This  will  affect  missile  performance  in  terms  of  time-to-go  and 
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higher  intercept  altitude.  A  quick  sensitivity  study  by  using  WIA  values  of  0,  100  and 
1000  and  computing  the  statistical  average  of  all  feasible  intercept  for  time-to-go, 
intercept  altitude  and  impact  angle  deviation  in  the  study  revealed  that  increasing  the 
weightage  by  a  factor  of  100  does  not  result  in  significant  changes  in  the  missile 
performance.  At  a  factor  of  1000,  the  effect  on  the  missile  flight  parameters  begins  to  be 
more  observable.  More  simulations  can  be  conducted  in  future  studies  for  higher  values 
of  the  WiA  that  will  cause  significant  effect  on  missile  performance.  Some  form  of  scaling 
factor  can  be  obtained,  from  which  a  missile  designer  can  select  suitable  values  Wia  to 
use  for  different  missions. 

There  is  hence  a  trade-off  between  the  weightage  coefficient  and  other 
performance  parameters  of  the  interceptor  missile.  From  the  trade-off  study,  it  was 
demonstrated  that  the  TS  guidance  allows  the  missile  designer  the  flexibility  to  adjust  the 
weightage  of  the  critical  parameters  in  order  to  optimize  the  intercept  path.  With  the 
insight  obtained  from  the  study,  the  effect  of  the  weightage  coefficient  on  missile 
performance  is  known.  There  is  a  need  to  balance  the  trade-off  between  the  weightage 
placed  on  the  different  critical  parameters.  In  the  particular  case  of  intercept  geometry, 
the  weightage  coefficient  (Wia)  must  be  carefully  selected  to  optimize  missile 
performance  for  a  particular  mission  or  scenario.  The  same  consideration  can  be 
translated  to  the  other  critical  parameters.  Again,  a  computer  program  can  be  developed 
to  aid  missile  designer  in  selecting  a  suitable  weightage  or  cost  function  coefficients  for 
different  mission  requirement  and  operational  scenario.  This  will  be  very  useful  if  TS 
guidance  in  implemented  in  intercept  missiles  and  offers  greater  advantage  over  other 
guidance  laws  as  one  that  allows  ‘customization’  for  different  mission  requirements 
simply  using  software  changes. 


D.  INTERCEPT  GEOMETRY  TRENDS 

The  simulation  study  also  provides  good  data  to  perform  simple  trending  analysis.  Useful 
charts  can  be  generated  to  provide  trending  information.  An  example  is  the  impact  angle 
deviation  data  derived  from  the  simulation.  It  provides  greater  insights  into  how  the 
intercept  geometry  will  change  with  respect  the  ballistic  missile  trajectory  azimuth  angle, 
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6bm,  with  respect  to  the  interceptor  launch  site.  The  angle  0bm,  is  defined  as  the  angle 
between  the  trajectory  plane  of  the  ballistic  missile  and  the  direction  from  the  interceptor 
launch  point  to  the  detected  ballistic  missile  position,  as  shown  in  Figure  26. 


Traiectorv 


A 


Ballistic 
Missile  Launch 
Direction 


Ballistic 
Missile 
Launch  Point 


Figure  26.  Illustration  of  Angle  between  interceptor  Launch  Position  and 

Ballistic  Missile  Trajectory  Plane 

Figure  27  shows  the  plot  of  impact  angle  deviation  for  interceptors  launched  from 
the  various  positions. 


Figure  27.  Effect  of  Interceptor  Launch  Direction  on  Intercept  Geometry 
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It  can  be  observed  from  the  plot  that  interceptors  launched  normal  to  the  ballistic 
missile  trajectory  plane  (0bm  =  90°)  usually  have  small  impact  angle  deviation  (very  close 
to  0°)  or  impact  the  ballistic  missile  at  close  to  90°  The  impact  angle  deviation  increases 
when  the  interceptor  flight  path  becomes  more  oblique  to  the  ballistic  missile  trajectory 
plane  (that  is,  Obm  becomes  smaller).  For  example,  the  impact  angle  deviation  is  larger, 
when  the  interceptor  is  fired  from  [80  km  N,  0  km  E],  as  compared  to  an  interceptor 
missile  being  launched  from  [0  km  N,  -50km  E] 

The  largest  impact  angle  deviation  occurs  when  the  interceptor  missile  is  in  the 
front-sector  of  the  ballistic  missile.  This  can  be  explained  from  the  fact  that  the  missile 
does  not  have  sufficient  time  to  effectively  ‘shape’  its  flight  path  for  a  90o  impact  on  the 
ballistic  missile,  compared  to  the  cases  where  the  interceptor  is  closing-in  from  the  side 
(normal  to  the  ballistic  missile  trajectory  plane).  By  plotting  the  results  for  the  other 
parameters  and  analyzing  them  in  future  studies,  other  useful  insight  into  the  interceptor 
missile  performance  can  be  deduced.  This  will  be  helpful  to  missile  designers  as  well  as 
operational  planners  to  optimize  the  potential  of  the  TS  guidance  for  possible  anti- 
ballistic  missile  defense  missions. 
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VI.  CONCLUSIONS 


The  simulation  study  verified  that  the  TS  guidanee  provides  feasible  solutions  to 
the  boost  phase  intercept  problem  involving  the  interception  of  a  fixed- trajectory  ballistic 
missile  by  a  surface-launched  interceptor  missile.  For  the  particular  scenario  used  in  the 
study,  it  was  shown  that  if  an  interceptor  system  using  the  TS  guidance  can  be  deployed 
within  30  to  60  km  (for  front  sector  intercept)  and  20  -  50  km  for  (rear-sector  intercept), 
it  is  possible  to  intercept  the  ballistic  missile  effectively  during  the  boost  phase.  The 
results  also  shows  that  the  new  guidance  allows  missile  designer  to  ‘customize’  the 
guidance  algorithm  for  specific  mission  by  changing  the  ‘weights’  of  some  critical 
parameters  However,  there  is  a  trade-off  between  the  each  of  the  weightage  coefficient 
and  missile  performance.  By  enhancing  the  weightage  on  intercept  geometry,  there  is  no 
change  in  the  time-to-intercept  or  intercept  altitude.  However,  the  region  for  the  feasible 
interceptor  launch  position  around  the  ballistic  missile  launch  point  has  reduced.  Hence, 
the  weightage  coefficient  will  have  to  be  carefully  chosen  to  ensure  that  there  is  a  balance 
of  benefit  and  penalty  with  respect  to  specific  missions.  This  aspect  of  ‘customization’ 
represents  a  clear  advantage  of  the  TS  guidance  over  other  guidance  law.  The  new 
guidance  algorithm  can  be  further  enhanced  and  fine-tuned.  Further  research  can  include 
conducting  simulation  study  using  a  6  DoF  model  of  the  TS  guidance,  combining 
guidance  solution  with  on-board  navigation  solution,  testing  the  complete  guidance, 
navigation  and  control  solutions  and  developing  a  computer  program  as  a  design  tool  to 
perform  trade-off  study  and  select  weightage  coefficients  for  different  mission 
requirements.  The  new  TS  guidance  holds  promise  to  be  a  feasible  and  implementable 
solution  to  the  boost  phase  ballistic  missile  intercept  problem  that  will  become  more 
prominent  in  the  coming  years. 
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APPENDIX  A:  TABLES  AND  PLOTS  OF  SIMULATION  RESULTS 


NO  OF  ITERATIONS 


WIA  =  0 


Ex  -  Exceed  Max  No.  of  Iterations  -  Non-feasible  Intercept 
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20 
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57 
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120 
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51.3 

47.1 

44.8 
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56.4 

50.3 

44.5 

40.5 

37.8 

80 
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53.1 

45.9 

39.5 

34.2 

30.8 

60 

100 

59 

50.7 

41.07 

34.5 

28.4 

23.9 

40 

100 

58 

48.2 

39.9 

31.3 

23.8 

17.6 

20 

100 

59 

48.3 

39.1 

30.3 

21.4 

0 

0 

100 

60 

50.7 

41 

31.8 

23.1 

0 

-20 

100 

100 

54.3 

45 

36.8 

28.8 

22.5 

-40 

100 
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61 

51.9 

43.5 

37 

32.3 

-60 

100 
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100 

60.2 

52.7 

47.6 

43.6 

-80 
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100 

59.2 

55.8 

Note  :  Value  of  0  and  100  Assigned  for  Non-Feasible  Intercept 
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IMPACT  ANGLE 


WIA  =  0 


Westing,  km 
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16 

18 
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27.3 
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80 

13 

15 

16.9 
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26.5 

28.9 

60 
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11 
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24.2 
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7.5 
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11.4 
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17.9 

26.5 

20 
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2.9 
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Ex 

0 

1.4 

3.2 

2.4 

2.5 

3 
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Ex 
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ALTITUDE 


WIA  =  0 


Westing,  km 
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100  or  Ex  -  Exceed  Max  No.  of  Iterations  -  Non-feasible  Intercept 
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NO  OF  ITERATIONS 


WIA=  100 


EX  -  No.  of  Iterations  Exceeded.  Represents  Non-Feasible  Region 
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IMPACT  ANGLE  DEF 
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INTERCEPT  ALTITUDE 
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IMPACT  ANGLE  DEVIATION 
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APPENDIX  B:  PARTICIPATION  IN  AIAA  YOUNG 
PROFESSIONAL  ROCKET  LAUNCH  COMPETITION  IN  MAY  2009 


INTRODUCTION 

The  author  participated  in  the  AIAA  Young  Professional  Rocket  Launch 
Competition  in  May  2009  as  a  member  of  the  NPS  Team  (see  Figure  28  for  the  team 
photo).  Team  Peacock,  one  of  the  two  teams  from  NPS,  comprised  eight  student 
members  and  a  faculty  advisor  (Prof  Yakimenko).  The  competition  required  its 
participant  to  apply  their  knowledge  to  build  a  rocket,  analyze  and  predict  it  performance 
and  launch  it  at  a  test  site,  using  common  rocket  kits,  a  choice  of  propulsion  units  and  a 
payload  provided  by  the  organizers.  The  8. 5 -feet  long  scaled  model  rocket  was  expected 
to  reach  height  in  excess  of  7,000  at  a  speed  of  greater  than  Mach  1 .  The  rocket  has  to  be 
recovered  after  reaching  its  maximum  height. 


Figure  28.  Team  Members  and  Advisor  with  Their  Rocket  at  the  Launch  Site 
RELEVANCE 

The  team’s  participation  in  the  competition  is  relevance  to  their  course  of  study  in 
NPS  in  general  in  that  it  provided  an  excellent  opportunity  for  the  students  to  apply  their 
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knowledge  in  aerodynamics,  propulsion,  missile  design  and  manufacturing  in  analyzing  a 
problems  and  physically  implementing  the  solutions  to  achieve  an  actual  task  or  a 
mission — in  this  case  building  the  rocket  and  launching  it.  In  the  process,  greater 
experience  and  knowledge  was  gained,  not  only  in  the  academic  areas,  but  also  in  terms 
of  project  management,  teamwork  and  communication.  This  is  a  good  example  of  the 
‘hands-on’  approach  towards  learning  and  an  interesting  and  impactful  means  of 
education.  With  respect  to  this  project,  this  rocket  is  like  a  scaled  model  of  an  interceptor 
missile  designed  to  intercept  a  ballistic  missiles  in  the  boost  phase.  The  characteristics  of 
the  rocket  resembles  that  of  a  hit-to-kill  interceptor,  vis  high  speed,  ‘light  weight’ 
(includes  only  the  body  and  fins,  actuators,  guidance  and  control),  small  payload  (perhaps 
a  small  KE  warhead  or  no  warhead  at  all)  and  ‘long  range’  (high  ratio  of  fuel  mass  to 
total  mass).  A  2-stage  rocket  would  bear  greater  resemblance  to  a  typical  missile 
interceptor.  The  valuable  experience  gained  provided  the  author  greater  insight  and  better 
understanding  of  the  fundamental  issues  facing  the  design  of  interceptor  missiles, 
particularly  the  aerodynamics,  propulsion,  stability  and  design  trade-offs. 


COMPETITION 

The  judging  criteria  for  the  competition  are  the  maximum  height  reached  by  the 
rocket  and  the  accuracy  of  the  predicted  performance  versus  actual  performance.  There  is 
a  winner  each  for  the  highest  height  reached  and  the  closest  predicted-to-actual 
performance  category.  Besides  the  actual  launch  itself,  a  report  submitted  by  each  team 
on  their  analysis  and  work  done  is  also  assessed  to  determine  the  winner  for  the  latter. 


GOALS 

Team  Peacock  set  itself  two  goals  for  this  competition  :  (1)  it  aimed  to  be  the 
winner  for  the  highest  altitude  category  and  (2)  to  build  the  most  stable  rocket  in  flight 
and  be  able  to  recover  the  rocket  successfully.  This  goal  guided  the  team  in  its  design  and 
construction  of  the  rocket. 
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TASKS 


The  team  took  about  4  months  to  build  the  rocket  and  prepare  for  the  launch, 
which  was  held  in  end  May  2009.  After  the  team  was  formed  in  Jan  09,  it  went  straight  to 
work  by  brainstorming  what  are  the  tasks  to  be  completed,  the  timeline,  resources 
required  (including  budget)  and  roles  of  individual  members.  Taking  into  account  the 
goals  and  tasks  identified,  allocation  of  tasks  and  a  schedule  broad  was  drafted. 


The  broad  tasks  facing  the  team  from  the  start  to  launch  were  as  follows: 

a.  Budgeting.  The  team  worked  within  a  budget  of  $1000  to  purchase  all  the 
required  materials  that  is  not  provided  by  the  organizer  (only  the  rocket  kit,  propulsion 
unit  and  the  common-use  payload  is  provided).  Table  8  provides  a  breakdown  of  the 
funds  needed  (for  supplies). 


Item 

Location 

Cost($) 

Qty 

Total($) 

60”  TAC-1  Parachute 

Giant  Leap  Rocketry 

86 

1 

86 

TAC-1  Kevlar  Bag 

Giant  Leap  Rocketry 

25 

1 

25 

Tubular  Kevlar  Shock  Cord  (1/2”  x  15') 

Giant  Leap  Rocketry 

25 

2 

50 

1/4”  Eyelet  Swivel 

Giant  Leap  Rocketry 

3 

4 

12 

AeroTech  Electronic  Fwd  Enclosure 

Aerotech-Rocketry 

170 

1 

170 

Construction  Materials 

Perfectflite.com 

100 

1 

100 

K270W(engine) 

Aerotech-Rocketry 

150 

2 

300 

Misc 

217 

TOTAL 


Table  8.  Budget  for  Rocket  Launch  Competition 


b.  Scheduling.  An  initial  schedule  was  drafted  (see  Figure  30)  so  that  members 
have  a  good  sense  of  the  key  milestones  and  timeline  for  the  various  preparatory 
activities.  The  schedule  was  updated  as  project  progresses  to  reflect  the  status  of  the 
various  tasks. 


59 


Taidlane 

Durdim 

RflcklLfflinth  Project 

IDifiys 

Simulition  f  Modeling 

Zlrf^s 

MatiiPteAdions 

21(!ays 

RoekiimPrtilicltoM 

Htntj  CticiMtti 

21tl>ys 

3itsvs 

Pufdukiffls 

Siisjfs 

14  dart 

AfflfflljIySMods 

Fdjys 

FdW 

RejglsAnolvsis 

FM  Report 

:(Jsyi 

Lauidi 

Idtrys 

Meeting 

jyjys 

Me 


Fefa15;QS  |F<!b22'09  Itol.tB  |MyC,TO 

s|m|t|w|Mf!s|s|m|t  |w|t!f  [sIsImIiIwItIf  |s|s|m|t|w|t|f|! 


IeImItIwIt 


Figure  29.  Timeline  for  the  Rocket  Launch  Project 


c.  Resource  Planning  and  Purchasing  (Supply).  Some  of  the  team  members 

were  assigned  the  role  to  plan  and  decide  the  items,  the  specifications  and 
quantities  needed  as  well  as  to  source  and  purchase  the  items  (the  actual 
purchasing  is  done  by  the  administrative  staff) 

d.  Structure  Design.  The  rocket  kit  is  provided  by  the  organizers.  It  uses  the 
Quasar  1/16  rocket,  which  comes  with  the  nose  cone,  body  sections,  fins  and  also 
a  parachute  for  recovery.  The  team  used  commercial  software  RockSim  to  predict 
the  CG  and  CP  calculation  as  well  as  to  determine  the  structural  dimensions  for 
the  rocket  design.  Figure  31  presents  the  CG  and  CP  position  and  the  static 
margin.  Based  on  the  calculation  and  prediction,  the  team  decided  to  reduce  the 
length  of  the  rocket  by  14  inches  as  well  as  halved  the  fin  area  to  reduce  the  static 
margin  to  2  and  optimize  the  performance  of  the  rocket  so  as  to  achieve  the 
required  stability  and  highest  altitude. 
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Quasar  Scale:  1/1  6 

Rocket  length:  101 .750  In. ,  diameter:  4.000  In. ,  span  diameter:  20.500  In. 
Rocket  mass  165.982  oz.  ,  Selected  stage  mass  165.982  gz. 

Shown  w/o  Engines. 


Method  CG  In.  CP  In.  CNa  Static  margin  Analysis 

Barrowman  G6.791  89.451  34.881  5.67  The  rocket  is  over  stable. 

RockSim  GG.791  90.307  42.1 43  5.88  The  rocket  is  over  stable. 

Figure  30.  The  Specifications  and  Dimensions  of  the  Quasar  Rocket  Kit 

e.  Rocket  Motor  Selection.  Two  engines,  both  with  about  the  same  impulse, 
were  provided  for  selection.  The  team  can  choose  either  the  engine  that  has  higher 
thrust  and  a  shorter  bum  time  (K800)  or  the  other  with  lower  thrast  but  longer 
bum  time  (K270).  After  conducting  simulation  study  of  the  predicted  flight 
profile,  based  on  the  rocket  modified  dimensions,  using  RockSim  and  a 
MATLAB  program  written  by  one  of  the  team  member,  it  was  found  that  the 
K270  best  meet  our  needs.  The  engine  motor  specifications  are  as  follows: 


Motor 

AeroTech  K270 

Diameter  (mm) 

54.0 

Length  (cm) 

57.9 

Prop.  Weight  (g) 

1,188.0 

Total  Weight  (g) 

2,100. 

Avg.  Thmst  (N) 

247.6 

Max.  Thmst  (N) 

425.7 

Tot.  Impulse  (Ns) 

2,154.9 

Burn  Time  (s) 

8.7 

The  motor  performance  chart  for  the  K270  are  shown  in  Figure  32. 
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Figure  3 1 .  The  Speeifications  and  Dimensions  of  the  Quasar  Roeket  Kit 

f.  Analysis  (Simulation  and  Modeling).  This  is  one  of  the  key  components  of 

this  rocket  launch  project.  The  team  has  to  use  various  methods  to  analyze  and 
predict  the  flight  performance/profile  of  the  rocket  to  be  built.  The  team  relied  on 
several  methods  to  analyze  the  flight  performance  of  the  rocket,  including  simple 
hand  calculations  for  rough  prediction,  using  RockSim,  a  commercial  rocket 
software  for  amateur  rocketeers  and  a  student-developed  MATLAB  program 
(with  Simulink  model)  written  by  one  of  our  members.  The  analysis  and 
simulation  done  helped  the  team  to  determine  the  structural  modification  needed, 
weight  and  balance,  CG  and  CP  position,  stability,  selection  of  engine  and  the 
maximum  altitude  prediction.  Figure  33  shows  the  Simulink  model.  Comparison 
of  results,  derived  from  RockSim  and  own  MATLAB  model,  can  be  made  and  it 
was  found  that  they  are  within  about  10%  of  each  other.  RockSim  predicted  a 
higher  maximum  height  reached  of  7,000  feet,  while  the  Simulink  model 
predicted  6,400  feet. 
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Figure  32.  The  Specifications  and  Dimensions  of  the  Quasar  Rocket  Kit 


g.  Testing.  Besides  analysis,  some  of  the  rocket  components  have  to  be 
physically  tested  to  have  a  higher  assurance  that  it  worked  as  designed  and  modifications 
made  did  not  introduce  problems  to  the  design.  The  team  tested  the  rocket  engine  at  the 
propulsion  laboratory  to  match  the  actual  performance  with  manufacturer  specifications, 
the  Electronic  Forward  Closure  (EFC)  to  test  the  trigger  mechanism  for  the  explosive 
charge  to  deploy  the  parachute  and  the  actual  parachute  to  select  the  one  to  be  used  for 
the  rocket  recovery  system.  During  the  first  propulsion  unit  test,  there  was  a  bum-through 
of  the  engine  casing  due  to  a  wrong  part  provided  by  the  manufacturer.  The  problem  was 
raised  to  the  manufacturer  and  a  replacement  part  was  delivered.  The  discovery  was  an 
important  one  as  it  prevented  the  same  occurrence  from  happening  during  the  actual 
launch  for  our  team  as  well  as  for  other  teams.  The  second  firing  was  very  successful  and 
it  verified  the  same  propulsion  characteristic  given  by  the  manufacturer.  The  EFC  test 
verified  that  the  trigger  mechanism  worked  and  would  cause  the  separation  of  the  rocket 
body  after  apogee  to  deploy  the  recovery  chute.  The  parachute  testing  helped  the  team  to 
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decide  using  a  single  chute  system  for  recovery  instead  of  a  dual-chute  system.  Figure  34 
shows  pictures  of  testing  being  conducted. 


Figure  33.  Team  Members  and  Advisor  with  Their  Rocket  at  the  Launch  Site 


h.  Construction.  The  actual  construction  of  the  rocket  took  about  three 

weeks,  after  all  the  parts  and  materials  were  delivered.  For  a  start,  a  rocket  stand  or  jig 
was  build  to  facilitate  the  assembling  work.  For  this  rocket,  the  body  was  made  up  of 
several  cylindrical  sections  of  thick  cardboard  material  and  the  nose  cone  is  made  of 
plastic.  The  construction  process  begins  with  strengthening  the  body  with  layers  of 
epoxy.  The  internal  compartments  and  fins  were  joined  to  the  body  using  fibre-glass  cloth 
with  epoxy  and  left  to  cure.  The  various  structures  were  attached  to  the  rocket  body  or 
inserted  inside  the  rocket  in  sections.  Fixing  the  fins  to  the  tail  section  required  more 
attention  to  ensure  they  are  symmetrically  attached.  Much  effort  was  put  in  to  ensure  the 
rigidity  as  the  fins  as  they  would  be  subjected  to  relatively  large  forces  in  flight.  The 
rocket  engine  casing,  EFC  and  parachute  were  then  put  inside  the  rocket  body.  Once  all 
the  components  were  assembled  in  their  sections,  the  three  sections  were  then  joined 
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together  using  internal  connectors  and  nose  cone  was  attached  to  the  body.  The  last  step 
is  painting  the  rocket  to  give  it  a  distinctive  color  for  aesthetic  reasons  as  well  as  to 
provide  a  smooth  finishing  to  reduce  skin  drag.  The  payload  would  be  placed  inside  the 
rocket  only  on  launch  day.  It  consisted  of  an  altimeter  that  would  send  altitude  data 
continuously  to  the  ground  receiver  for  recording  purposes.  Figure  37  shows  some 
pictures  of  the  rocket  being  built  at  the  laboratory. 


Figure  34.  Team  Members  and  Advisor  with  Their  Rocket  at  the  Launch  Site 


i.  Launch.  The  launch  was  the  climax  of  the  five-month  effort  and  the  team 
was  as  prepared  as  it  could  be.  It  was  held  on  26  May  09,  Saturday  at  Koehn  Lake 
Launch  Site  near  to  Mojave.  On  that  day,  there  were  9  teams  present  with  their  rockets. 
Team  Peacock  was  the  first  team  to  launch.  After  the  rocket  propellant  was  inserted  into 
the  casing,  it  was  mounted  onto  the  launch  rail  for  the  ignition.  The  igniter  was  fired  by 
the  firer,  who  was  positioned  at  the  firing  ‘bunker’.  All  other  observers  were  herded  to 
the  observation  deck  to  watch  the  firing.  At  the  countdown  of  10,  the  rocket  was 
launched.  It  lifted-off  vertically  with  a  sharp  ‘bang’  and  within  seconds,  the  rocket  motor 
accelerated  the  rocket  to  more  than  3,000  ft  and  out  of  visual  sight,  leaving  a  trail  of 
smoke  behind  it.  Moments  later,  the  parachute  was  deployed  and  the  rocket  descended  to 
the  ground  in  two  parts,  as  expected.  It  landed  ‘safely’  at  the  desert  ground  close  to  the 
launch  site  and  was  retrieved  by  the  recovery  team.  Figure  36  shows  the  launch  of  the 
rocket. 
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Figure  35.  Team  Members  and  Advisor  with  Their  Rocket  at  the  Launch  Site 


RESULTS 

The  results  made  all  the  hard  work  pay  off.  Team  Peacock  emerged  the  winner  of 
the  closest  altitude  prediction  category.  The  recorded  maximum  altitude  reached  during 
the  actual  launch  was  5,700  feet.  This  was  about  9%  deviation  from  the  predicted  height 
based  on  the  team’s  analysis.  Figure  37  presents  the  results  in  terms  of  the  altitude  versus 
time  plot  that  was  recorded  by  the  organizer. 
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Time,  s 


Figure  36.  Official  Altitude  Plot  for  Team  Peacock 


LESSON  LEARNED 

The  end  result  is  not  as  important  as  the  process  of  getting  it.  Though  the  team  did 
not  expect  to  win  and  entered  the  competition  with  an  altitude  to  experience  and 
understand  rocket  design  and  rocketry  better,  the  results  was  nonetheless  a  good  one  for 
first-timers  like  us.  There  were  many  lessons  learnt  from  the  success  and  failure  of  some 
teams.  By  analyzing  and  discussing  the  causes  of  failure  with  the  other  teams,  it  provided 
much  insight  into  the  considerations  that  went  into  the  design  and  some  of  the  cardinal 
errors  that  were  made.  Using  the  example  of  a  rocket  that  ‘exploded’  mid-air,  their  over- 
ambitious  fm  design,  which  reduced  too  much  of  the  control  surface  for  normal  flight  and 
stability,  it  provided  a  good  lesson  for  our  team  as  well  in  our  future  design.  In  future 
participation,  more  prediction  tools  and  program  should  be  made  available  for  team 
members  to  analysis  the  aerodynamics  and  flight  performance.  The  project  may  be  part  of 
the  NPS  missile  design  course,  culminating  in  the  actual  launch  of  the  rocket;  all  these 
within  2  quarters.  It  can  help  to  reinforce  some  of  the  fundamental  principles  of 
aerodynamics  and  rocket  design.  In  summary,  this  was  a  worthwhile  project  that 
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combined  learning,  application  of  knowledge  and  great  fun  into  one,  and  certainly  helped 
in  providing  a  better  grasp  of  fundamentals  in  missile  design. 
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APPENDIX  C:  INDUCED  DRAG  PROGRAMME  DEVELOPED 


1.  ZLDragC_mod.m 


function  [CD]  =  ZLDragC_mod (M, load_f , mass, rho, vel, Sref ) 

%  Written  by  LTC  Weng  Wai  Leong,  Naval  Postgraduate  School,  Dec  2009 

%  This  function  interpolates  the  known  drag  polars  for  various  Mach 
number  to  determine  drag  coefficient.  Input  variables  are  Mach  number, 

%  load  factor,  mass,  density,  velocity,  referenece  area  the  interceptor 
or  ballistic  missile  (BM)  to  be  utilized 

%  in  the  BMFlightS.m,  SMTrajectory  and  SMFlightS.m  programs  (written  by 
LT  Lukacs)  for  the  calculation  of  forces  acting  on  the 
interceptor/Ballistic  Missile. 

%  The  CL  and  Cd  data  are  input  in  Excel  files.  Data  obtained  from  Prof 
Yakimenko 

%  Variable  List 

%  M  =  mach  number  of  interceptor/BM 
%  load_f=  load  factor  in  the  normal  plane 
%  mass  =  mass  of  interceptor/BM 
%  rho  =  density  of  air 
%  vel  =  velocity  of  interceptor/BM 
%  Sref  =  reference  area  of  interceptor/BM 

global  CL_data  Cd_data 

Cd_curve  =  interpl (Cd_data ( : , 1 ) , Cd_data ( : , 2 : 11 ) , M,  ' nearest '  ,  ' extrap ' ) ; 
%interpolate  the  Cd  data  in  the  array  for  the  new  Mach  number 

CL_curve  =  interpl (CL_data ( : , 1 ) , CL_data ( : , 2 : 11 ) , M, ' nearest ' , ' extrap ' ) ; 
%interpolate  the  CL  data  in  the  array  for  the  new  Mach  number 

p  =  polyf it (CL_curve, Cd_curve, 2 ) ;  %  generate  a  new  drag  polar  (CL_Cd 

curve)  for  new  Mach  number  using  2nd  order  polynominal  function 

k=p(l)  ; 

CL0=0; 

CD0=p (3) -k^CL0^2; 

CL= (2^1oad_f ^mass^9 . 81) / ( rho^vel^2 ^Sref ) ; 

if  abs(CL)>5,  CL^S'^sign  (CL)  ;  end 

CD  =  CDO  +  k^ (CL-CLO) ^2; 
return 
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APPENDIX  D:  MODELING  PROGRAM  CODE 


A.  3DOF  INTERCEPTOR 
1.  SMFlightS.m 


%%  This  is  the  Main  M-File  for  the  simulation.  Run  this  file  with  all 
the  required  function  files  in  the  same  current  directory. 

%  Written  by  LT  John  A.  Lukacs  IV,  Naval  Postgraduate  School,  June  2006 
%  Corrected  by  O.Yakimenko,  August  2007,  November  2009 

%  This  script  develops  and  tracks  the  flight  path  of  the  interceptor 
%  missile.  For  the  first  ten  seconds  it  integrates  a  series  of 
%  acceleration  commands  to  simulate  a  vertical  launch  and  tip  over. 
Upon 

%  activation  of  the  guidance  law,  it  sends  the  known  values  to  the 
guidance 

%  law  and  receives  back  the  future  time  history  of  the  optimal  flight 
path . 

%  This  script  then  implements  that  optimal  path.  It  updates  the  final 
%  conditions  and  recalculates  the  optimal  flight  path  at  an  interval  of 
10 

%  seconds.  The  script  calls  BRFlight3,  SMParamsS.m,  SMDrag . m, STatmos . m, 

%  and  SMGuidance.m 


%%  List  o 
%  Acc_SM 
%  AllSM 
%  alt 
%  CD 
%  count 
interval 
%  dist 
%  Drag 
%  Forces 
%  g 


f  variables 

=  index-based  vector  of  acceleration  values 
=  index  based  vector  of  all  interceptor  values 
=  altitude 
=  drag  coefficient 

=  counting  variable  to  determine  guidance  law 


update 


SM 


=  cumulative  distance  travelled 
=  total  drag  force 

=  index-based  vector  of  force  values 
=  gravitational  force,  based 


on  WGS-84  value  of 


gravitational 

o 

o 

%  m_i  = 

%  Model_SM 
%  MV 
sound) 

%  N1,N2 
length 

%  num_SM  = 

%  nx,ny,nz 
respectively 
%  nx (y, z ) _Op  = 
%  path  = 

%  Pos_SM 
%  press  = 


attraction  and  altitude 
interceptor  mass 

index-based  vector  of  internal  values 

=  interceptor  speed  in  Mach  (relative  to  local  speed  of 

=  variables  used  to  ensure  optimal  vectors  are  the  same 

number  of  interations  conducted  (used  for  plotting) 

=  axial  acceleration  command,  body  frame  x,  y  and  z, 

index  based  vector  of  the  optimal  flight  path  values 
returned  time  history  of  the  optimal  path 
index-based  vector  of  position  values 
local  atmospheric  pressure 
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%  psi 
%  psidot 
%  psidot_Op 
%  px,py,pz 
%  px (y, z)_Op 
%  py_Op 
%  pz_Op 

%  q 
%  Re 
%  ro 
%  Sref 
%  state 
%  t 

%  target 
%  temp 
%  tgo 
%  th 
%  th_Op 
%  thdot 
%  thdot_Op 
%  Thrust 
%  time_Op 
%  time_SM 
%  update 
%  V 

%  V_Op 
%  Vdot 
%  Vdot_Op 
%  Vel_SM 
%  w 

interceptor 


heading  angle 

rate  of  change  of  heading  angle 

index  based  vector  of  the  optimal  flight  path  values 

X,  y  and  z  components  of  position 

index  based  vector  of  the  optimal  flight  path  values 

index  based  vector  of  the  optimal  flight  path  values 

index  based  vector  of  the  optimal  flight  path  values 

counting  variable  (used  in  another  program) 

WGS-84  Earth's  radius 
local  atmospheric  density 

planar  reference  area  (for  drag  calculations) 
the  state  of  the  interceptor 
global  time 

the  state  of  the  rocket 
local  atmospheric  temperature 
time  to  go  to  intercept 
flight  path  angle 

index  based  vector  of  the  optimal  flight  path  values 
rate  of  change  of  flight  path  angle 

index  based  vector  of  the  optimal  flight  path  values 
thrust  generated  by  interceptor  motor 

index  based  vector  of  the  optimal  flight  path  values 
index-based  vector  of  time  values 
number  of  updates  to  the  guidance  law  conducted 
velocity  of  the  interceptor 

index  based  vector  of  the  optimal  flight  path  values 
rate  of  change  of  the  velocity 

index  based  vector  of  the  optimal  flight  path  values 
index-based  vector  of  velocity  values 

=  number  of  0.5s  steps  between  the  lunches  of  traget  and 


close  all,  clear  all,  clc 


global  Re  alphag  path  %  shared  by  SMFlightS,  SMGuidance  & 

SMTra j  ectory 

global  q  states  %  shared  by  SMFlight3  &  SMGuidance 

global  Pos_BR  Pos_SM  Npol  Npt  kpolar  weight90  %  ...  by  SMFlightS  & 

SMTra j  ectory 

%%  Computing  the  flight  path  of  a  Ballistic  Missile  in  a  gravity  turn 
BRFlightS 

%%  Initializing  variables  for  an  Interceptor 
t=0;  dt=0.5;  q=0;  w=120;  %  60  sec  detection  time 
Nupd=30;  count=Nupd;  update=0; 

Npt=100;  %  the  number  of  points  the  optimal  trajectory  is 

computed  at 

Npol=8;  %  number  of  coefficients  in  approximating 

polynomials 

kpolar=0 ; 

prompt  =  { 'Northing  wrt  to  BM  launch  point,  km' , . . . 

'Westing  wrt  to  BM  launch  point,  km', . . . 
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'Weighting  coefficient  for  impact  angle'}; 
dlg_title  =  'Enter  Interceptor  launch  coordinates'; 
num_lines  =  1; 
def~=  {  '  80  '  ,  '  60  '  ,  '  1000  '  }  ; 

answer  =  inputdlg (prompt , dlg_title, num_lines , def, ' on ') ; 
px  =str2num  (answer  { 1 })  "^1000  ; 

py  =str2num(answer{2})'^1000; 

weight 90=str2num (answer { 3 } ) ; 

%  px=-100000/20; 

%  py=100000/2; 

%  weight 90=100 ; 
pz=Re ; 

px_old=px;  py_old=py;  pz_old=pz; 
dist=0 ; 

psi=atan2 ( Pos_BR (120,2) -py, Pos_BR (120, 1 ) -px) ;  psi_old=psi ; 
th^OO'^pi/lSO ;  th_old=th; 

V=l;  V_old=V; 

%%  Computing  Interceptor  flight  path 
for  i=l:20%1000 
t=t+dt ; 

%%  Boost  phase  (vertical  launch) ,  t<10s 
if  t<10 

%  Speed,  Mach  number 
alt=norm ( [px; py ; pz ] ) -Re  ; 

[ ro, press , temp] =STatmos (alt) ; 

MV  =  V/sqrt (1 . 402^287 . 053^temp) ; 

%  Forces 

g=3 . 98 6004 4 18e 14 /norm ( [px;py;pz] ) ^2; 

[m_i,Sref]  =  SMParams3 (t) ; 

CDtable  =  SMDrag (MV) ; 

if  t<=6.  Thrust  =  206000;  psidot=0; 
else  Thrust  =  95300;  psidot=0; 

end 

ny  =  V/g'^cos  ( th)  "^psidot ; 
nz  =  V/g'^thdot+cos  ( th)  ; 
nn  =  sqrt (ny^2+nz^2 ) ; 

CL=  ( 2 '^nn'^m_i'^g)  /  ( ro'^V^2 "^Sref )  ; 

CD  =  CDtable (1)  +  kpolar^CL^2; 

Drag  =  ro^V^2 ^CD^Sref /2  ; 
nx= (Thrust-Drag) /m_i/ g; 

alphag=180/pi'^nn'^m_i'^g/  (ro'^V^2'^Sref /2 )  /13; 

%  Kinematics 
Vdot^g"^  (nx-sin  (th)  )  ; 
psidot^ny^g/V/cos (th) ; 
thdot=g^ (nz-cos (th) ) /V; 

%  Collecting  variables 
t ime_SM ( i , 1 ) =t ; 


thdot=0 ; 

thdot=- 


0.075; 
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Forces_SM (i, 1) =nx; 
Model_SM(i, 1) =V; 


Forces_SM (i, 2 ) =ny; 
Model_SM ( i , 2 ) =Vdot ; 
Model_SM ( i , 4 ) =thdot ; 
Model_SM (i, 6) =psidot  ; 
Pos_SM (i, 2 ) =py; 


Forces_SM (i, 3) =nz; 


Model_SM(i, 3) =th; 
Model_SM (i, 5) =psi; 
Pos_SM (i, 1) =px; 


Pos_SM (i, 3) =pz ; 


Pos_SM(i, 4) =dist; 

Vel_SM (i, 1) =V^cos (th) ^cos (psi) ; 

Vel_SM  (i,  2 )  =V'^cos  (th)  "^sin  (psi)  ; 

Vel~SM(i, 3) =V^sin (th) ; 

Acc_SM (i, 1) =Vdot^cos (th) ^cos (psi) -V^cos ( th) ^sin (psi) ^psidot . . . 

-V'^sin  (th)  "^cos  (psi)  "^thdot; 

Acc_SM  (i,  2 )  ^Vdof^cos  (th)  "^sin  (psi)  +V'^cos  (th)  "^cos  (psi)  "^psidot .  .  . 

+V'^sin  (th)  "^sin  (psi)  "^thdot; 

Acc_SM  (i,  3)  ^Vdot'^sin  (th)  +V'^cos  (th)  "^thdot ; 

%  Euler  integration 
V= V_o 1 d+ Vdo  t^dt; 
psi=psi_old+psidot^dt  ; 
th=th_old+thdot^dt  ; 
px=px_old+V'^cos  (th)  "^cos  (psi)  "^dt; 
py=py_old+V^cos (th) ^sin (psi) ^dt; 
pz=pz_old+V'^sin  (th)  "^dt; 

dist= (dist+abs (norm ( [px-px_old; py-py_old; pz-pz_old] ) ) ) ; 

V_old=V;  psi_old=psi;  th_old=th; 

px_old=px;  py_old=py;  pz_old=pz; 

else 

%%  Optimal  guidance,  t>10s 

if  count==Nupd  &  update==0  %  Recomputing  the  trajectory  every 

Nupd  cycles 


update=update+l ; 

fprintf (' Starting  Interceptor '' s  Guidance  Update 


#%2 . Of \n ' , update) 


state= [px;py;pz; V; th;psi; Vdot; thdot ; psidot ] ; 
target^ [ Pos_BR ( i+w, 1 ) ; Pos_BR ( i+w, 2 ) ; Pos_BR ( i+w, 3 ) ; 


Vel_BR(i+w, 1) ;Vel_BR(i+w, 2) ; Vel_BR ( i+w, 3 ) ; 
Acc_BR ( i+w, 1 ) ; Acc_BR ( i+w, 2 ) ; Acc_BR ( i+w, 3 ) ] ; 


else 


N2=length (path  (:,!)); 


end 

%  Identifying  variables 

time_Op ( : , update) = [path ( : , 1 ) ; zeros (N1-N2 , 1 ) ] ; 
px_Op ( : , update) = [path ( : , 2 ) ; zeros (N1-N2 , 1 ) ] ; 
py_Op ( : , update) = [path ( : , 3 ) ; zeros (N1-N2 , 1 ) ] ; 
pz_Op ( : , update) = [path ( : , 4 ) ; zeros (N1-N2 , 1 ) ] ; 

V_Op ( : , update) = [path ( : , 5 ) ; zeros (N1-N2 , 1 ) ] ; 
th_Op ( : , update) = [path ( : , 6) ; zeros (N1-N2 , 1 ) ] ; 
psi_Op ( : , update) = [path ( : , 7 ) ; zeros (N1-N2 , 1 ) ] ; 
Vdot_Op ( : , update) = [path ( : , 8 ) ; zeros (N1-N2 , 1 ) ] ; 
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thdot_Op ( : r update) = [path ( : , 9) ; zeros (N1-N2 , 1 ) ] ; 
psidot_Op ( : , update) = [path ( : , 10 ) ; zeros (N1-N2 , 1 ) ] ; 
nx_Op ( : , update) = [path ( : , 8 ) ; zeros (N1-N2 , 1 ) ] ; 
ny_Op ( : , update) = [path ( : , 9) ; zeros (N1-N2 , 1 ) ] ; 
nz_Op ( : , update) = [path ( : , 10 ) ; zeros (N1-N2 , 1 ) ] ; 
count=l ; 

end 

if  count==101  %  added  by  OY  2009-10-08 

count=count-l ; 
elseif  count==0 
count=l ; 

end 


t=time_Op (count , update) ; 
nx=nx_Op (count, update)  ; 
ny=ny_Op (count, update)  ; 
nz=nz_Op (count , update) ; 

V=V_Op (count , update) ; 

Vdot=Vdot_Op (count , update) ; 
th=th_Op (count, update) ; 

thdot=thdot_Op (count , update) ; 
psi=psi_Op (count , update)  ; 

psidot=psidot_Op (count, update) ; 
px=px_Op (count, update) ; 
py=py_Op (count, update) ; 
pz=pz_Op (count , update) ; 


Forces_SM (1,2) =ny  ; 
Model_SM (1,2) =Vdot ; 
Model_SM (1,4) =thdot ; 
Model_SM (1,6) =psldot ; 
Pos_SM (1,2) =py; 


%  Collecting  variables 
t lme_SM ( 1 , 1 ) =t ; 

Forces_SM (1,1) =nx; 

Model_SM(l, 1) =V; 

Model_SM(l,  3)  =th; 

Model_SM ( 1 , 5 ) =psl ; 

Pos_SM (1,1) =px; 

Pos_SM(l, 4) =dlst; 

Vel_SM (1,1) =V^cos (th) ^cos (psl)  ; 

Vel_SM  (1,2)  =\['^cos  (th)  "^sln  (psl)  ; 

Vel_SM(l, 3) =V^sln (th) ; 

Acc_SM (1,1) =Vdot^cos (th) ^cos (psl) -V^cos (th) ^ 
-V^sln (th) ^cos (psl) ^thdot; 

Acc_SM (1,2) =Vdot^cos ( th) ^sln (psl) +V^cos (th) ^ 
+V'^sln  (th)  "^sln  (psl)  '*'thdot; 

Acc_SM  (1,3)  ^Vdot'^sln  (th)  +V'^cos  (th)  "^thdot ; 
Update (1,1) ^update; 


Forces_SM ( 1 , 3 ) =nz ; 


Pos_SM (1, 3) =pz ; 


sin  (psl) ^psldot .  .  . 
cos (psl) ^psldot. . . 


end 


%  Time  step 
count=count+l ; 

dlst= (dlst+abs (norm ( [px-px_old; py-py_old; pz-pz_old] ) ) ) ; 

V_old=V; 

psl_old=psl ; 

th_old=th; 

px_old=px; 

py_old=py; 

pz_old=pz ; 

%  the  end  of  the  "If"  loop 
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end  %  the  end  of  the  "for"  loop 

A11SM=  [time_SM  Forces_SM  Model_SM  Pos_SM  Vel_SM  Acc_SM] ;  %  update]  ; 


2.  SMParam3.m 


function  [SM  mass. 

Sref ] =SMParams3 (t) 

o 

o 

Written  by  LT 

John  A.  Lukacs  IV,  Naval  Postgraduate  School 

,  June  2006 

o 

o 

This  function 

calculates  the  J  Matrix  and  Mass  of  the  intercepter. 

o 

o 

assuming  a  cruciform  rocket  in  two  stages  to  intercept.  This 

o 

o 

function  also 

returns  the  reference  (base)  diameter  of  the 

missile . 

o 

o 

%  Notes: 

o 

o 

The  first  stage 

last  6  seconds,  the  second  stage  lasts  an 

additional 

o 

o 

10  seconds. 

Stage  1  (booster)  seperates  upon  completion. 

Stage  2 

o 

o 

does  not  separate  after  completion 

o 

O  ■ 

%  List  of  variables 

o 

o 

dia 

reference  diameter,  base  diameter 

o 

o 

1 

length,  varies  by  component 

o 

o 

p  SM  stl  fuel 

density  of  stage  1  rocket  fuel 

o 

o 

p  SM  st2  fuel 

density  of  stage  2  rocket  fuel 

o 

o 

p  SMstr 

density  of  structural  material 

o 

o 

r 

radius,  varies  by  component 

o 

o 

ro 

outer  radius,  varies  by  component 

o 

o 

ri 

inner  radius,  varies  by  component 

o 

o 

SM  mass 

total  rocket  mass 

o 

o 

SM  nose 

total  mass  of  nosecone  section 

o 

o 

SM  stl  fcr 

consumption  rate  of  stage  1  fuel 

o 

o 

SM  stl  fuel 

remaining  stage  1  fuel  based  on  time  and 

o 

o 

consumption  rate 

o 

o 

SM  stl  str 

total  mass  of  stage  1  structural  material 

o 

o 

SM_stl_tfm 

total  mass  of  stage  1  fuel 

o 

o 

SM_st2_fcr 

consumption  rate  of  stage  2  fuel 

o 

o 

SM  st2  fuel 

remaining  stage  2  fuel  based  on  time  and 

o 

o 

consumption  rate 

o 

o 

SM_st2_str 

total  mass  of  stage  2  structural  material 

o 

o 

SM_st2_tfm 

total  mass  of  stage  2  fuel 

o 

o 

t 

time 

o 

o 

th 

structural  thickness 

o 

o 

V  bo  dy 

=  volume  of  body  structural  material 

o 

o 

V  nose  str 

volume  of  nosecone  structural  material 

o 

o 

V  nose  strO 

volume  of  nosecone  structural  material. 

o 

o 

intermediate  value 

o 

o 

V  nose  strl 

volume  of  nosecone  structural  material. 

o 

o 

intermediate  value 

o 

o 

V  stl  fuel 

volume  of  stage  1  fuel 

o 

o 

V  stl  str 

volume  of  stage  1  structural  material 
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%  V_st2_fuel  =  volume  of  stage  2  fuel 

%  V_st2_str  =  volume  of  stage  2  structural  material 

%%  Structural  components 
p_SMstr  =  4225; 
th  =  .0208; 

%  Mass  of  nose  cone 
1  =  .8255; 
r  =  0.34/2; 

V_nose_str0  =  pi^ ( 1^ ( ( r^2+1^2 ) / ( 2 ^r ) ) ^2-1^3/3- ( ( ( r^ 
(2^r) ) -r) ^ ( (r ^2  +  1^2) / (2^r) ) ^2^asin  (1/  (  (r^2  +  l^ 
1  =  .8255-th; 
r  =  0.34/2-th; 

V_nose_strl  =  pi^ ( 1^ ( ( r^2+1^2 ) / ( 2 ^r ) ) ^2-1^3/3- ( ( ( r^ 
(2^r) ) -r) ^ ( (r ^2  +  1^2) / (2^r) ) ^2^asin (1/  (  (r^2  +  l^ 
V_nose_str  =  V_nose_str0-V_nose_strl+pi'^r^2  "^th; 
SM_nose  =  1 . 3^V_nose_str^p_SMstr ; 

%  Mass  of  Body/Warhead  Section 
1  =  .849; 
ro  =  0.34/2; 
ri  =  0.34/2-th; 

V_body  =  I'^pi'^  ( ro^2-ri^2  )  ; 

SM_body  =  V_body'^p_SMstr+115  ; 

%  Mass  of  Stage  1  (Mk72  Booster) 

1  =  1.72; 
ro  =  0.53/2; 
ri  =  0.53/2-th; 

V_stl_str  =  l'^pi'^(ro^2-ri^2)+2'^pi'^ro^2'^th; 
SM_stl_str  =  V_stl_str'^p_SMstr  ; 

%  Mass  of  Stage  2  (Mkl04  Engine) 

1  =  2.88-2^th; 
ro  =  0.34/2; 
ri  =  0.34/2-th; 

V_st2_str  =  I'^pi'^  ( ro^2-ri^2 ) +2 '^pi'^ro^2 "^th; 
SM_st2_str  =  V_st2_str'^p_SMstr  ; 

%%  Fuel  Components 

%  Stage  1  Solid  Fuel  is  HTPB/AP/Al 

1  =  1.72; 

ri  =  0.53/2-th; 

V_stl_fuel  =  0 . 80'^l'^pi'^ri^2  ; 
p_SM_stl_f uel  =  1860; 

SM_stl_tfm  =  468; 

SM_stl_fcr  =  468/6; 

%  Stage  2  Solid  Fuel  is  TP-H1205/6 

1  =  2.88; 

ri  =  0.34/2-th; 

V_st2_fuel  =  0 . 60'^l'^pi'^ri^2  ; 

SM_st2_tfm  =  360; 

p_SM_st2_f uel  =  SM_st2_tfm/V_st2_f uel ; 


2+1^2) / . . . 
2)/(2^r)))); 


2  +  1^2) /  .  .  . 
2)/(2^r)))); 
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SM_st2_fcr  =  360/15; 
if  t<6 

%%  Stage  1  -  Stage  1  Fuel  is  consumed.  Stage  2  Fuel  is  not  used. 
SM_stl_fuel  =  SM_stl_tfm  -  SM_stl_fcr  ^  t; 

SM_mass  =  SM_nose+SM_body+SM_st2_str+SM_st2_tfm+SM_stl_str+ . . . 

SM_stl_fuel ; 
dia  =  0.53; 

elseif  t<21 

%%  Stage  2  -  Stage  1  has  separated.  Stage  2  Fuel  is  consumed. 
SM_st2_fuel  =  SM_st2_tfm  -  SM_st2_fcr  ^  (t-6) ; 

SM_mass  =  SM_nose+SM_body+SM_st2_str+SM_st2_f uel ; 
dia  =  0.34; 

else 

%%  Stage  3  -  The  unpowered  nosecone  and  Stage  2  remains . 

SM_mass  =  SM_nose+SM_body+SM_st2_str ; 
dia  =  0.34; 

end 

Sref=pi'^dia^2/4  ; 
return 


3.  SMDrag.m 


Written  by  LT  John  A.  Lukacs  IV,  Naval  Postgraduate  School,  June  2006 
and  renamed  by  Prof  Oleg  Yakimenko  in  November  2009 

%  This  function  interpolates  to  determine  drag  coefficient  based  on  the 
%  Mach  number  and  the  boost  or  glide  phase  of  the  rocket  to  be  utilized 
%  in  the  BMFlight.m  and  SMFlight.m  programs  for  the  calculation  of 
%  forces  acting  on  the  rocket/missile. 

%  The  tables  used  were  point-plotted  from  Prof  Hutchins'  ME4703 
%  "Missile  Flight  Analysis"  Class  Notes 

%%  Setting  a  List  of  Variables 

%  BDrag  =  boost  phase  drag  interpolation  table 
%  GDrag  =  glide  phase  drag  interpolation  table 
%  M  =  mach  number 

%  Mach  =  mach  number  interpolation  table 
%%  Tables: 

Mach  =  [0  0.90  1.1  1.2  1.5  2.0  2.5  3.0... 

3.5  5.0  6.0]; 

BDrag  =  [0.1444  0.1444  0.2778  0.2778  0.2308  0.1778  0.1481  0.1296... 
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0.1185  0.1000  0.0950] ; 

GDrag  =  [0.2461  0.2461  0.4615  0.4615  0.3615  0.2846  0.2500  0.2192... 

0.2000  0.1500  0.1300]  ; 

%plot (Mach, BDrag, ' o-- ' , Mach, GDrag, ' +- . ' ) 

%%  Calculation  of  Drag  Coefficient  Values: 
if  M  >  6 
M  =  6; 

end 

BDrag  =  interpl (Mach, BDrag, M, ' cubic ') ; 

GDrag  =  interpl (Mach, GDrag, M, ' cubic ') ; 

Drag= [BDrag; GDrag] ; 
return 


B.  3DOF  TARGET 


1.  BRFUght3.m 


o 

o 

%  Complementary  M-File 

o 

o 

Written 

by 

LT  John  A.  Lukacs  IV,  Naval  Postgraduate 

School,  June  2006 

o 

o 

This  script  integrates  the  position. 

velocity. 

and  acceleration 

values 

o 

o 

at  each 

time  step  to  determine  the  flight 

path  of 

a 

ballistic  missile 

o 

o 

in  a  gravity  turn.  The  script  calls  BRParams3.m, 

ZLDragC.m,  and 

o 

o 

STatmos . 

m 

o  ^ 
O  ■ 

%  Setting  a  : 

List  of  Variables 

o 

o 

acc 

total  acceleration 

o 

o 

alt 

altitude 

o 

o 

ax 

X  component  of  acceleration 

o 

o 

ay 

y  component  of  acceleration 

o 

o 

az 

z  component  of  acceleration 

o 

o 

CD 

drag  coefficient 

o 

o 

Drag 

total  drag  force 

o 

o 

dt 

time  step  interval 

o 

o 

g 

=  gravitational  force. 

based 

on 

WGS-84  value  of 

gravitational 

o 

o 

attraction  and  altitude 

o 

o 

gm 

initial  launch  angle 

o 

o 

i 

interval  count 

o 

o 

m  r 

rocket  mass 

o 

o 

Mspd 

rocket  speed  in  Mach  (relative 

to  local 

speed  of  sound) 

o 

o 

num  BR 

number  of  interations  conducted  (used 

for  plotting) 

o 

o 

nx 

axial  force 

o 

o 

press 

local  atmospheric  pressure 

o 

o 

px 

X  component  of  position 

o 

o 

PY 

y  component  of  position 

o 

o 

pz 

z  component  of  position 
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%  Re 
%  ro 
%  spd 
%  Sref 
%  t 

%  temp 
%  Thrust 
%  vx 
%  vy 
%  vz 

%  Acc_BR 
%  Forces_BR 
%  Pos_BR 
%  time_BR 
%  Vel_BR 
%  All  BR 


WGS-84  value  for  Earth's  radius 
local  atmospheric  density 
rocket  speed  in  m/s 

planar  reference  area  (for  drag  calculations) 
time 

local  atmospheric  temperature 
thrust  generated  by  motor 
X  component  of  velocity 
y  component  of  velocity 
z  component  of  velocity 

index-based  vector  of  acceleration  values 
index-based  vector  of  force  values 
index-based  vector  of  position  values 
index-based  vector  of  time  values 
index-based  vector  of  velocity  values 
index  based  vector  of  all  rocket  values 


%%  Initializing  Variables 
dt=0.5;  i=0; 
Re=6.378137e6; 


%%  Setting  initial 
px=0 ; 
py=0; 
pz=Re ; 

gm=75^pi/180 ; 
vx=cos (gm) ; 
vy=0; 

vz=sin (gm) ; 


Conditions 


%%  Computing  Ballistic  Flight  Path 
for  t=0:dt:200%22%19.5 
i=i+l; 


%  Speed,  Mach  Number 
spd=norm( [vx;vy;vz] ) ; 
alt=norm ( [px; py ; pz ] ) -Re ; 
if  alt<86000 

[ ro, press , temp] =STatmos (alt )  ;  %  Calling  STatmos  function 

else 

[ ro, press , temp] =STatmos ( 8 6000 ) ;  %  Calling  STatmos  function 

end 

Mspd=spd/ sqrt  ( 1 . 4 02 "^2 87 temp)  ; 


%  Forces 

g=3 . 98 6004 4 18e 14 /norm ( [px;py;pz] )^2; 

[m_r , Sref ] =BRParams3 ( t ) ;  %  Calling  BRParamsS 

function 

CD=BRDrag (Mspd) ;  %  Calling  BRMrag  function 


if  t<125 

Thrust^lOSOOO"^ 9 . 81  ; 
CD=CD (1) ; 
elseif  t<240 
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Thrust=29950^9.81; 

CD=CD (1) ; 

else 

Thrust=0 ; 

CD=CD (2) ; 

end 

Drag=ro^spd^2^CD^Sref /2  ; 
nx= (Thrust-Drag) /m_r/g; 

%  Accelerations 

g=3 . 98  6004  418el4'^pz/norm  (  [px;py;pz]  )  ^3; 
ax^g'^nx'^cos  (gm)  ; 
ay=0  ; 

az^g'^nx'^sin  (gm)  -g; 
acc=norm( [ax;ay;az] ) ; 

%  Collect  Variables 
time_BR ( 1 , 1 ) =t ; 

Pos_BR (1,1) =px; 

Pos_BR (1,2) =py; 

Pos_BR ( 1 , 3 ) =pz ; 

Pos_BR (1,4) =norm ( [px; py ; pz ] ) ; 

Vel_BR (1,1) =vx; 

Vel_BR (1,2) =vy ; 

Vel_BR (1,3) =vz ; 

Vel~BR(i, 4) =spd; 

Vel_BR(i,  5)  =Mspd; 

Acc_BR (1,1) =ax ; 

Acc_BR (1,2) =ay; 

Acc_BR ( 1 , 3 ) =az  ; 

Acc_BR (1,4) =acc ; 

Acc_BR (1,5) =nx  ; 

Forces_BR (1,1) =Thrust; 

Forces_BR (1,2) =m_r ; 

Forces_BR (1,3) =Drag; 

%  Time  Step 
px^px+df^vx; 
py=py+dt^vy; 
pz^pz+df^vz  ; 
vx^vx+df^ax; 
vy^vy+df^ay ; 
vz^vz+df^az ; 
num_BR=length (time_BR) ; 

end 

%plot (time_BR( : , 1) , (Pos_BR( : ,3) -Re) /lOOO) 

A11BR=  [time_BR  Forces_BR  Pos_BR  Vel_BR  Acc_BR] ; 

clear  CD  Drag  Mspd  Sref  Thrust  acc  alt  ax  ay  az  dt 

clear  g  gm  h  1  m_r  nx  press  px  py  pz  ro  spd  t  temp  vx  vy  vz 
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2.  BRParams3.m 


function  [BR_mass, Sref , length] =BRParams (t) 

%  Written  by  LT  John  A.  Lukacs  IV,  Naval  Postgraduate  School,  June  2006 

%  This  function  calculates  the  mass  of  the  rocket,  assuming  a  cruciform 
%  rocket  in  two  stages  plus  an  unpowered  nosecone  stage.  This  function 
%  also  returns  the  reference  (base)  diameter  of  the  missile. 


%%  Notes: 

%  The  first  stage  last  125  seconds,  the  second  stage  lasts  an 
additional 

%  110  seconds.  The  stages  seperate  upon  completion. 


%%  Setting  a  List  of  Variables 
%  BR_mass  =  total  rocket  mass 
%  BR_nose  =  total  mass  of  nosecone  section 
%  BR_stl_bt  =  burntime  for  stage  1 

%  BR_stl_fcr  =  consumption  rate  of  stage  1  fuel 

%  BR_stl_fuel  =  remaining  stage  1  fuel  based  on  time  and  consumption 


rate 

%  BR_stl_str 
%  BR_stl_tfm 
%  BR_st2_bt  = 
%  BR_st2_fcr 
%  BR_st2_fuel 
rate 

%  BR_st2_str 
%  BR_st2_tfm 
%  dia 


=  total  mass  of  stage  1  structural  material 

=  total  mass  of  stage  1  fuel 

burntime  for  stage  2 

=  consumption  rate  of  stage  2  fuel 

=  remaining  stage  2  fuel  based  on  time  and  consumption 

=  total  mass  of  stage  2  structural  material 

=  total  mass  of  stage  2  fuel 

=  reference  diameter,  base  diameter 


%  t  =  time 


%%  Setting  Structural  Components 
BR_nose  =  250; 

BR~st2_str  =  2288; 

BR_stl_str  =  9000; 

%%  Setting  Fuel  Components 
BR_stl_tfm  =  50970; 

BR~stl~bt  =  125; 

BR_stl_fcr  =  BR_stl_tfm/BR_stl_bt; 

BR_st2_tfm  =  12912; 

BR_st2_bt  =  110; 

BR_st2_fcr  =  BR_st2_tfm/BR_st2_bt; 
if  t<125 

%%  Stage  1  -  Stage  1  Fuel  is  consumed.  Stage  2  Fuel  is  not  used 
BR_stl_fuel  =  BR_stl_tfm-BR_stl_f  cr'^t  ; 

BR_mass  =  BR_nose+BR_stl_str+BR_stl_f uel+BR_st2_str+BR_st2_tfm; 
dia  =  2.2; 
length  =  2+14+16; 

elseif  t<240; 

%%  Stage  2  -  Stage  1  has  separated.  Stage  2  Fuel  is  consumed 
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BR_st2_fuel  =  BR_st2_tfm-BR_st2_fcr^ (t-125) ; 

BR_mass  =  BR_nose+BR_st2_str+BR_st2_f uel ; 
dia  =  1.3; 
length  =  2+14; 

else 

%%  Stage  3  -  Stage  2  has  separated,  only  the  unpowered  nosecone  remains 
BR_mass  =  BR_nose; 
dia  =  1.3; 
length  =  2; 

end 

Sref =pi'^dia^2 /  4  ; 
return 


3.  BRDrag.m 


function  Drag  =  BRDrag (M) 

%  Written  by  LT  John  A.  Lukacs  IV,  Naval  Postgraduate  School,  June  2006 

%  Renamed  by  Oleg  Yakimenko,  November  2009  to  be  used  as  separate  drag 
file 

%  for  interceptor  and  ballistic  missile  instead  of  the  common  ZLDragC.m 
file 

%  used  previously 

%  This  function  interpolates  to  determine  drag  coefficient  based  on  the 
%  Mach  number  and  the  boost  or  glide  phase  of  the  rocket  to  be  utilized 
%  in  the  BMFlight.m  and  SMFlight.m  programs  for  the  calculation  of 
%  forces  acting  on  the  rocket/missile. 

%  The  tables  used  were  point-plotted  from  Prof  Hutchins'  ME4703 
%  "Missile  Flight  Analysis"  Class  Notes 

%%  List  of  variables 


%  BDrag  =  boost  phase  drag  interpolation  table 
%  GDrag  =  glide  phase  drag  interpolation  table 
%  M  =  mach  number 

%  Mach  =  mach  number  interpolation  table 


%%  Tables: 

Mach  =  [0  0.90 

3.5  5.0 

BDrag  =  [0.1444  0.1444 
0.1185  0.1000 
GDrag  =  [0.2461  0.2461 
0.2000  0.1500 


1.1  1.2  1.5 

6.0]; 

0.2778  0.2778  0.2308 
0.0950] ; 

0.4615  0.4615  0.3615 
0.1300] ; 


.0  2.5  3.0. . . 

.1778  0.1481  0.1296. . . 
.2846  0.2500  0.2192. . . 
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%plot (Mach, BDrag, ' o-- ' , Mach, GDrag, ' . ' ) 

%%  Calculation  of  drag  coefficient  for  the  boost  and  glide  phases: 
if  M  >  6 
M  =  6; 

end 

BDrag  =  interpl (Mach, BDrag, M, ' cubic ') ; 

GDrag  =  interpl (Mach, GDrag, M, ' cubic ') ; 

Drag= [BDrag; GDrag] ; 
return 


C.  GUIDANCE  ALGORITHMS 
I.  SMGuidance.m 


function  path=SMGuidance (time, state, target, i) 

%  Written  by  LT  John  A.  Lukacs  iV,  Naval  Postgraduate  School,  June  2006 
%  Corrected  by  O.Yakimenko,  August  2007,  October  2009 

%  This  function  takes  in  the  state  of  the  interceptor  and  target  and 
%  generates  an  initial  guess  at  the  final  conditions  (position, 

%  orientation  angles,  range,  and  time  to  intercept)  through  a  first- 
order 

%  trajectory  assumption  and  iterative  process.  it  then  calls  the 
%  fminsearch  function  using  those  initial  guesses.  Finally,  it  plots 
the 

%  returned  optimal  flight  path  and  associated  variables. 


%%  List  o 
%  best 

o 

o 

%  BC 
%  cost 
%  costs 
iteration 
%  free 

o 

o 

%  init 
%  J 
%  N 

%  nmax 

o 

o 

%  path 
%  psi 
%  psidot 
%  psif 

o 

o 

%  psit 
%  Py,Pz 
%  q 


f  variables 

=  vector  of  the  variables  in  the  optimal  path  returned  from 
the  fminsearch  function 
=  boundary  conditions 

=  cost  function  value  returned  from  SMGuidanceCost  function 
=  array  of  the  value  of  the  cost  variables  at  each 

=  variables  that  fminsearch  can  modify,  specifically 
[  tau; thf ; psif] 

=  vector  of  initial  estimates 
=  vector  of  cost  function  variable  values 
=  length  of  the  path  vector  (used  for  plotting) 

=  maximum  acceleration  capability  of  the  interceptor, 
altitude  dependent 

=  returned  time  history  of  the  optimal  path 
=  initial  interceptor  heading  angle 

=  initial  rate  of  change  of  interceptor  heading  angle 
=  final  interceptor  heading  angle,  calculated  from  final 
conditions  estimate 
=  target  heading  angle 

=  penalty  functions  on  the  y  and  z  acceleration 
=  variable,  counting  trajectory  updates  during  intercept 
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%  qq 

optimization 
%  range  = 

%  state  = 

o 

o 

%  states  = 

%  target  = 

o 

o 

%  tau_f  = 

%  tgo  = 

%  th 

%  thdot  = 

%  thf 

%  tht  = 

%  tic . . toe  = 

%  trys  = 

%  V 
%  V_f 
%  Vave 
%  Vdot 
%  xO 

%  xdO  = 

%  xddO 
%  xdf 

%  xdt  = 

%  xf  = 

%  xmult  = 

%  xt  = 

%  ymult  = 

%  zmult  = 


=  variable,  counting  the  number  of  iterations  during 

estimate  of  distance  between  target  and  interceptor 

state  of  the  interceptor  missile,  sent  from  SMGuidance.m 

[px; py; pz ; V; th; psi ; Vdot ; thetadot ; psidot ] ; 

array  of  the  values  of  all  processes  in  SMGuidance.m 

state  of  the  rocket,  sent  from  SMGuidance.m  at  i+w, 

synchronizing  times 

value  of  the  virtual  arc 

time  to  go  to  intercept 

initial  interceptor  flight  path  angle 

initial  rate  of  change  of  interceptor  flight  path  angle 

final  interceptor  flight  path  angle 

target  flight  path  angle 

MATLAB  function  to  track  run  time 

vector  of  optimal  path  and  derivative  values 

initial  interceptor  velocity 

final  interceptor  velocity 

average  interceptor  velocity 

initial  interceptor  acceleration 

initial  inteceptor  position 

initial  interceptor  velocity 

initial  inteceptor  acceleration 

final  inteceptor  position 

current  target  velocity 

final  inteceptor  acceleration 

ratio  value  (used  for  plotting) 

current  target  position 

ratio  value  (used  for  plotting) 

ratio  value  (used  for  plotting) 


global  Re  alphag  path  %  shared  by  SMFlightS,  SMGuidance  & 

SMTra j  ectory 

global  q  states  %  shared  by  SMFlightS  &  SMGuidance 

global  qq  thdot  psidot  tgo  costs  trys  %  ...  by  SMGuidance  & 

SMTra j  ectory 


%%  Counting  trajectory  updates  during  intercept 

q=q+l; 


%%  Initializing  Interceptor  (all  states) 
xO=state (1:3)  ; 

V=state (4) ; 
th=state ( 5 ) ; 
psi=state ( 6 ) ; 

Vdot=state (7) ; 
thdot=state (8)  ; 
psidot=state (9)  ; 
xdO  =  [V^cos (th) ^cos (psi) ; 

V'^cos  (th)  "^sin  (psi)  ; 

V^sin ( th) ] ; 

xddO  =  [Vdof^cos  ( th)  "^cos  (psi) -V'^cos  ( th)  "^sin  (psi)  "^psidot- 

V^sin  (th) ^cos (psi) ^thdot; 

Vdot'^cos  (th)  "^sin  (psi)  iV'^cos  (th)  "^cos  (psi)  "^psidot- 
V^sin  (th)  "^sin  (psi)  thdot; 
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Vdot'^sin  (th)  +V'^cos  (th)  "^thdot]  ; 

%%  Initializing  BM  (up  to  the  second-order  derivatives) 
xt  ^target (1:3) ; 
xdt  ^target ( 4 : 6 ) ; 
xddt=target (7:9)  ; 

%%  Estimating  time-to-go 
tgol=100;  delta=5; 
vmaxhyp=2  60  0; 
while  delta>l 
if  tgol>20 

V_f=vmaxhyp-10'^  (tgol-2  0)  ; 

Vave=  (2  0'^vmaxhyp/2+  (tgol-20)  ( vmaxhyp+V_f )  /2)  /tgol; 

else 

V_f=tgol'^vmaxhyp/2  0  ; 

Vave=V_f/2; 

end 

xf ^xt+xdf^tgol ; 

tgo2=sqrt (  (xf  ( 1 ) -xO  ( 1 ) ) "2+ (xf ( 2 ) -xO ( 2 ) ) "2+ (xf ( 3 ) -xO ( 3 ) ) "2  )  .  .  . 

/ (norm (xdt) +Vave) ; 
delta=abs (tgo2-tgol) ; 
tgol= ( tgol+tgo2 ) /2  ; 

end 

tgo=tgol ; 

%%  Initializing  optimization 

range=sqrt ((xf(l)-x0(l))"2+((xf(2)-x0(2)))"2+((xf(3)-x0(3)))"2); 

%fprintf ( ' Tra j ectory  update  #  %2.0f  \n',q) 

fprintf ( ' \nSlant  range  to  target:  %5.1fkm  \n ' , range/10^3 ) 

tau_f  =0 . 00045'^range-1000'^  (q-1 )  ;  %  guess  on  tau_f 

fprintf (' Guess  on  the  virtual  arc  length:  %5.1f  \n',tau_f) 
fprintf (' Guess  on  the  time-to-go:  %5.1fkm  \n',tgo) 

predicted_xdt=xddt'^tgo ; 

tht  =atan2 (predicted_xdt ( 3 ) , norm (predicted_xdt (1:2) ) ) ; 
psit=atan2 (predicted_xdt (2) , predicted_xdt (1) ) ; 
thf  =0;%-tht;  %  guess  on  thf 

psif  =psi ; %psit+pi ;  %  guess  on  psif 

free  = [tau_f; thf ; psif ; . 1 ; 1 ; 1 ; -0 . 001 ; -0 . 001 ] ; 

BC= [xO ; xdO ; xddO ; time ; xt ; xdt ; xddt ] ; 

%%  Searching  for  the  minimum  performance  index 
qq=0;  %  counting  iterations  to  converge 

tic 

options=optimset ( ' Maxi ter ',100, 'Tolfun',1, 'TolX',1) ; 

best  =  fminsearch ( @ (x)  SMTraj ectory (x, BC) , free, options ) ;  % 

Optimization 

tcpu=toc; 

fprintf (' \nlt  took  %6.0f  interations  to  converge\n ' , qq) 

fprintf (' Elapsed  time  is  %6.1f  seconds\n ' , tcpu) , 

fprintf (' Combined  performance  index  is  %5 . If \n ',  costs (end,  1 ) ) 

fprintf ( '  including:  tau_f=%5.1f,  t2go=%5.1fs,  ImpAngle=%4 . If 

off,\n',... 

costs  (end, 2 ) , costs (end, 3) , 180 /pi^ costs (end,  4 )  ) 
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fprintf ( '  Pny=%5.1f  and 

Pnz=%5.1f\n' , costs (end, 5) , costs (end,  6)  ) 

%fprintf ( '  Dtgo=%5.1f  and 

Altf ine=%5 .lf\n' , costs (end, 7 ) , costs (end, 8 ) ) 

tau_f=best ( 1 ) ; 
thf  =best (2 ) ; 
psif  =best ( 3 ) ; 

[best (4) ; best (5) ; best (6) ; best (7) ;best (8) ] ; 

V_f=path (end,  5)  ; 
xf =path (end, 2:4) ; 

fprintf ([' Impact  occurs  at  Altitude  of  %4.1fkm,  Northing=%4 . If km,  and 

I 

' Westing=%4 . If km\n ' ] ,  [xf ( 3 ) -Re  xf ( 1 ) 

xf (2) ] 710^3) 

%fprintf  (  '  with  interceptor '' s  speed  of  %4.1f  km/s\n  '  ,  V_f /10"^3 ) 

xdf=  [V_f "^cos  (thf)  "^cos  (psif)  ; 

V_f "^cos  (thf)  "^sin  (psif)  ; 

V_f^sin (thf)  ]  ; 

%%  Plotting  results 

%{ 

xmult=20000/ (norm (xdO ) +norm (xdt ) ) ; 
ymult=20000/ (norm (xdO ) +norm (xdt) ) ; 
zmult=20000/ (norm (xdO ) +norm (xdt) ) ; 
nmax=40+ (40-10) / (0-50000) ^ (path ( : , 4) -Re) ; 

f igure (' Name Bird-eye  view')  %  Bird-eye  view 

plot3  (path ( : , 2) 710^3, path ( : , 3) 710^3,  (path ( : , 4) -Re) 710^3,'- 
. b ' , ' Linewidth ' , 2 ) 
hold  on,  grid 

plot3  (10'"-3^  [xO  (1)  -xmult^xdO  (1)  ;x0  ( 1 )  +xmult^xd0  (1)  ]  ,  .  .  . 

10^-3"^  [xO  (2)  -ymulf^xdO  (2)  ;x0  (2)  +ymult'^xd0  (2)  ]  ,  .  .  . 

10^-3^ [xO (3) -Re-zmult^xdO (3) ;x0 (3) - 
Re+zmulf^xdO ( 3 ) ] , ' c ' , ' Linewidth ' , 2 ) 

plot3  (xf  (1)  7l0'"3,xf  (2)  710'"3,  (xf  (3)  -Re)  7l0'"3,  'pr  '  ,  '  Linewidth ',  2 ) 
plot3  (10^-3"^  [xf  (1)  -xmulf^xdt  (1)  ;xf  (1)  +xmult'^xdt  (1)  ]  ,  .  .  . 

10^-3"^  [xf  (2)  -ymulf^xdt  (2)  ;xf  (2)  +ymult'^xdt  (2)  ]  ,  .  .  . 

10^-3"^  [xf  (3)  -Re-zmulf^xdt  (3)  ;xf  (3)  - 
Re+zmulf^xdt ( 3 ) ] , ' r ' , ' Linewidth ' , 2 ) 

plot3  (xO  (1)  7l0'"3,x0  (2)  710^3,  (xO  (3)  -Re)  710^3,  '  ^b'  ,  '  Linewidth ',  5 ) 
plot3  (  (xO  (1)  +xmult^xd0  (1)  )  7l0'"3,  (xO  ( 2 )  +xmult^xd0  (2)  )  710^3,  .  .  . 

(x0(3)- 

Re+xmulf^xdO  (3))710"^3,  '"^c',  '  linewidth  '  ,  2 ) 

plot3 ( (xf (1) +xmult^xdt (1) ) 710^3, (xf ( 2 ) +xmult^xdt (2) ) 710^3, . . . 

(xf(3)- 

Re+xmulf^xdt  (3))710'^3,  '^r',  '  linewidth  '  ,  2 ) 

hl=legend (' Intercept  tra j ectory Interceptor '' s  velocity  vector',... 

'Impact  point ',' BM '' s  velocity  vector  at 

intercept ' , ' Location ' , ' Best ' ) ; 
set  (hi,  ' Fonts ize ' , 8 ) ; 

xlabel (' Northing  (km)'),  ylabel (' Westing  (km)'),  zlabel (' Altitude 
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(km)  ') 

%title (' Interception  Geometery', ' Fontsize ' , 10 ) 
view (-102,8) 
axis  equal 

f igure (' Name Combined  PI  and  Virtual  arc')  %  Performance  index  & 

Virtual  arc 
subplot  (211 ) 

semi logy (costs(:,l)/costs(l,l) , 'h-. ') ,  grid 

xlabel (' Iteration ') ,  ylabel (' Relative  PI  (PI_i/PI_l') 

xlim([l  qq] ) ,  axis  'auto  y'  %  ylim([0  2]) 

subplot (212 ) 

plot  (costs (:, 2 ),' h- .')  ,  grid 

xlabel  (' Iteration ') ,  ylabel (' Length  of  virtual  arc,  \it\tau_f') 
xlim( [1  qq] ) 

f igure (' Name ',' Delta  Time-to-go  and  Altitude  violation')  %  Dtgo  & 
Altitude  fine 
subplot  (211 ) 

plot (costs (:, 7 ),' h- .') ,  grid 
xlim ( [1  qq] ) 

xlabel (' Iteration ')  ,  ylabel (' \Delta  \itt_{go}  \rm(s)') 
subplot (212) 

plot  (costs (:, 8) /10^3,  ' h- .  '  )  ,  grid 

xlabel  (' Iteration ') ,  ylabel ('Alt.  violation,  (km)') 
xlim( [1  qq] ) 

f igure  (' Name ',' Impact  angle  and  Time-to-go')  %  Impact  angle  &  Time- 

to-go 

subplot (211 ) 

plot (real ( 180/pi^acos  (costs ( : , 4 ) ) ) ,  ' h- .  ' ) ,  grid 
axis ( [ 1  qq  60  90 ] ) 

xlabel  (' Iteration ')  ,  ylabel  (' Impact  angle  ('^o)') 
subplot (212) 

plot  ( costs (:, 3 ),' h- .')  ,  grid 

xlabel (' Iteration ') ,  ylabel (' Time-to-go,  \itt_{go}  \rm(s)') 
xlim( [1  qq] ) 

f igure (' Name ',' G-load  factors')  %  G-load  constraints 

subplot (211 ) 

plot  (path ( : , 1 ) , path ( : , 9 ) ,  ' -b .  ' ) ,  grid 
hold  on 

plot  (path ( : , 1 ) , path (:,10) ,  '--g.  ') 
plot (path ( : , 1 ) , nmax ( : ) , ' r ' , ' Linewidth ' , 2 ) 
plot (path ( : , 1 ) , -nmax ( : ) , ' r ' , ' Linewidth ' , 2 ) 
hl=legend ( ' n_y ' , ' n_z ' , ' Dynamic  constraints ' , 2 ) ; 
set  (hi ,  ' Font Size ' , 8 ) ; 

xlabel ('Time  (s)'),  ylabel ('Load  factor  (g) ' ) 
subplot (212) 

plot(costs(:,5)/10^7, '-b. ', ' Linewidth ' , 2 ) ,  grid 
hold  on 

plot  (costs  (:,  6)  /  lO'^l ,  '  --g .  '  ,  '  Linewidth  '  ,  2 ) 
xlabel  (' Iteration ') ,  ylabel (' Relative  penalty') 
hl=legend ( ' n_y  penalty ' , ' n_z  penalty ' , ' Location ' , ' Best ' ) ; 
set  (hi,  ' Fonts ize ' , 8 ) ; 


88 


xlim ( [1  qq] ) 

f igure (' Name Lambda  and  Tau  profile')  %  Lambda  and  tau 

subplot (211 ) 

plot  (path (i  fl)  f path ( : , end) ,'.'),  grid 
xlabelCTime  (s)'),  ylabel  (  ' \it\lambda  '  ) 
subplot (212 ) 

plot  (path (:  fl)  f path ( : , end-1 ) ,'.'),  grid 
xlabel ( ' Time  (s) ' ) ,  ylabel ( ' \it\tau ' ) 

f igure (' Name Speed  and  Angle  of  attack  profile')  %  SM  speed  and  Angle 
of  attack 
subplot  (211 ) 

plot (path (i  fl)  f path (:, 5 ),'.') ;  grid 
xlabel ('Time  (s)'),  ylabel (' Speed,  V  (m/s)') 
subplot (212 ) 

plot  (path ( : , 1) , alphag,  '  .  ' ) ,  grid 

xlabel  ('Time  (s)'),  ylabel  (' Angle  of  attack  ("^o)') 

f igure (' Name ',' Euler  angles  profile')  %  SM  Euler  angles 

subplot (211 ) 

plot  (path ( : , 1 ) , path (:, 6)^180 /pi,  '  .  ',  ' Linewidth ' , 2 ) ,  grid 
hold  on 

plot  (path  (end,  1)  ,thf'^180/pi,  'ro') 
ylim( [-30  90] ) 

xlabel ('Time  (s)'),  ylabel (' \theta  (^o)') 
subplot  (212 ) 

plot (path ( : , 1 ) , path (:, 7)^180 /pi, 'g. ', ' Linewidth ' , 2 ) ,  grid 
hold  on 

plot  (path  (end,  1)  ,psif'^180/pi,  'ro') 
ylim( [-180  180] ) 

XlabelCTime  (s)'),  ylabel  (' \psi  ('"o)') 

%} 

%%  Creating  results  structure 

states { q, 1 } =path; 

states { q, 2 } =BC; 

states {q, 3 }=free; 

states {q, 4 }=best; 

states { q, 5 } =costs ; 

return 


2.  SMGuidanceCost.m 


function  [cost, J, Py, Pz] =SMGuidanceCost (free, const) 

%  Written  by  LT  John  A.  Lukacs  IV,  Naval  Postgraduate  School,  June  2006 

%  This  function  calculates  the  cost  of  the  proposed  trajectory  returned 
%  from  the  SMTra j ectory . m  function  based  on  the  optimization  parameters 
%  and  penalty  parameters  defined  herein.  This  is  a  sub-funtion  of  the 
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%  SMGuidance . m  function's  fminsearch .  This  cost  value  is  used  to 
%  determine  whether  the  proposed  trajectory  is  optimal.  The  trajectory 
%  that  returns  the  minimum  value  of  J  is  the  optimal  function. 


%%  Variable 
%  calccost 

o 

o 

%  const 

o 

o 

%  cost 
function 
%  costs 
%  dist 
%  free 

o 

o 

%  init 
%  J 
%  N 

%  nmax 

o 

o 

%  nx 
%  ny 
%  nz 
%  path 

o 

o 

nz  '  ] 

%  psi 
%  psidot 
%  psif 

o 

o 

%  Py 
%  Pz 

%  qq 
%  t 

%  tau_f 
%  tgo; 

%  th 
%  thdot 
%  thf 
%  time 
%  V 
%  V_f 
%  Vdot 
%  X 
%  xO 
%  xdO 
%  xddO 
%  xdf 
%  xdt 
%  xt 


List 

=  a  global  variable  of  the  value  of  the  J  function  (used 
for  plotting) 

=  variables  that  fminsearch  cannot  modify,  including  system 
contraints,  specifically  [xO ; xdO ; xddO ; time ; xt ; xdt ; init ] 

=  cost  function  value  returned  from  SMGuidanceCost . m 

=  vector  of  values  of  the  cost  variables  at  each  iteration 
=  cumulative  distance  travelled 

=  variables  that  fminsearch  can  modify,  specifically 
[  tau; tgo; thf; psif] 

=  vector  of  initial  estimates 
=  vector  of  cost  function  variable  values 
=  length  of  the  path  vector  (used  for  plotting) 

=  maximum  acceleration  capability  of  the  interceptor, 
altitude  dependent 

=  axial  acceleration  command,  body  frame  x 
=  axial  acceleration  command,  body  frame  y 
=  axial  acceleration  command,  body  frame  z 
=  returned  time  history  of  the  optimal  path,  specifically 
[time'  X(l:3, :) '  V  th '  psi'  Vdot'  thdot'  psidot'  nx '  ny' 

=  initial  interceptor  heading  angle 

=  initial  rate  of  change  of  interceptor  heading  angle 
=  final  interceptor  heading  angle,  calculated  from  final 
conditions  estimate 

=  penalty  function  on  the  y  acceleration 

=  penalty  function  on  the  z  acceleration 

=  counting  variable 

=  current  time 

=  value  of  the  virtual  arc 

=  time  to  go  to  intercept 

=  initial  interceptor  flight  path  angle 

=  initial  rate  of  change  of  interceptor  flight  path  angle 
=  final  interceptor  flight  path  angle 
=  optimal  path  time  history 
=  initial  interceptor  velocity 
=  final  interceptor  velocity 
=  initial  interceptor  acceleration 

=  the  optimal  path  time  history  in  cartesian  coordinates 
=  initial  inteceptor  position 
=  initial  interceptor  velocity 
=  initial  inteceptor  acceleration 
=  final  inteceptor  velocity 
=  current  target  velocity 
=  current  target  position 


[path] =SMTra j ectory (free, const)  ; 


global  calccost  costs  q  qq  tgo  trys 
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%  Initialize  Variables 

qq=qq+l; 

dist=0 ; 

Re=6.378137e6; 

%%  Identify  Variables 
time=path ( : , 1 ) ; 

X=path ( : , 2 : 4 ) ; 

V=path  (  : , 5 )  ; 
th=path  (  : , 6 ) ; 
psi=path  (  \  fl) ; 

Vdot=path  (  : , 8 )  ; 
thdot=path ( : , 9) ; 
psidot=path ( : , 10) ; 
nx=path  (  : , 11 )  ; 
ny=path  (  : , 12 )  ; 
nz=path  (  : , 13 ) ; 

N=length (path  (  : , 1) ) ; 

tau_f=f ree  ( 1 ) ; 
tgo=f ree  (2 )  ; 
thf=f ree  ( 3 ) ; 
psif=f ree  ( 4 ) ; 

x0=const  (1:3) ; 
xd0=const (4:6) ; 
xdd0=const  (7 : 9) ; 
t=const  ( 10 )  ; 
xt=const  (11 : 13) ; 
xdt=const (14:16) ; 
init=const (17 : 19) ; 

V_f=path (N, 5) ; 

xdf= [V_f ^cos (thf ) ^cos (psif ) ; 

V_f '*'cos  (thf)  "^sin  (psif)  ; 

V~f^sin (thf)  ]  ; 
tgo=path (N, 1 ) ; 

for  i=2 : 1 : N 

dist=dist+abs (norm( [X(i,l)-X(i-l,l) ;X(i,2)-X(i-l,2)  ;X(i,3)-X(i- 
1.3);])); 
end 

nmax=40+ (40-10) / (0-50000) ^ (path ( : , 4) -Re) ; 

%%  Calulate  Cost  of  the  chosen  trajectory 
J= [  tau_f; 
tgo; 

lOO'^abs  (dot  (xdf ,  xdt )  / norm  (xdf )  / norm  (xdt )  )  ]  ; 

Py=sum (max ( 0 , abs (ny) -nmax) .^2) ; 

Pz=sum (max ( 0 , abs (nz ) -nmax) .^2) ; 

cost=0 . 3 3 ones  (1,3)  J+norm  (  [  Py ;  Pz  ]  )  ; 
costs (qq, 1:6)= [cost; J; Py; Pz ; ] ; 
calccost=cost ; 
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return 


3.  SMTrajectory.m 


function  cost=SMTra j ectory (free, BC) 

%  This  function  computes  a  candidate  trajectory  and  associated  cost 
based  on 

%  the  vector  of  varied  parameters  "free"  and  boundary  conditions  "BC" . 

%  This  is  a  sub-funtion  of  the  SMGuidance.m  function's  fminsearch. 

%  This  function  creates  a  7th  order  set  of  equations  and  evaluates  that 
set  at 

%  the  boundary  conditions  supplied  by  the  inputs.  It  then  calculates 
the  time 

%  history  of  all  the  flight  vehicle  variables,  including  controls  and 
reactions , 

%  necessary  to  develop  that  flight  path.  A  plot  command  set  at  the  end 
of  this 

%  function  will  plot  a  chart  of  the  iterations  at  the  end  of  run  if 
desired . 

%  Finally,  this  function  calculates  the  cost  of  the  candidate 
tra j  ectory 

%  combining  the  value  of  the  performance  index  and  penalties.  This  cost 
value  is 

%  used  to  determine  whether  the  proposed  trajectory  is  optimal.  (The 
tra j  ectory 

%  that  returns  the  minimum  value  of  the  cost  is  the  sub-optimal  one.) 

%  O.Yakimenko,  Naval  Postgraduate  School,  November  2009 


%%  List  of  variables 

%  A, Ax, Axp, Axpp, Axppp  =  cell  matrices  of  coefficients  of  a  candidate 
reference 

%  trajectory  and  their  derivatives  wrt  virtual 


arc 

%  BC  =  boundary  conditions,  specifically 

[xO ; xdO ; xddO ; time ; xt ; xdt ; xddt ] 

%  Cx, Cxp, Cxpp, Cxppp  =  coefficients  of  a  candidate  reference 

trajectory  and 

their  derivatives  wrt  virtual  arc 
=  tau  step  value 
=  time  step  value 

=  variable  parameters,  specifically  [ tau; thf ; psif ] 

=  gravitational  force 
=  lambda,  virtual  speed 

=  first-order  derivative  of  lambda  wrt  to  virtual  arc 
=  maximum  acceleration  capability  of  the  interceptor, 
altitude  dependent 

%  nX, nXp, nXpp, nXppp  =  norm  of  reference  trajectory  and  its 


%  dtau 
%  dtime 
%  free 

%  g 
%  L 
%  Lp 
%  nmax 

o 

o 
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derivatives 

%  nx  =  axial  acceleration  command,  body  frame  x 

%  ny  =  axial  acceleration  command,  body  frame  y 

%  nz  =  axial  acceleration  command,  body  frame  z 

%  path  =  returned  time  history  of  the  optimal  path,  specifically 

%  [time'  X(l:3,:)'  V  th '  psi'  nx '  ny'  nz '  tau '  L'] 

%  psi  =  interceptor  heading  angle 

%  psidot  =  rate  of  change  of  interceptor  heading  angle 

%  qq  =  variable  counting  the  number  of  iterations 

%  t  =  global  current  time 

%  tau  =  virtual  arc 

%  tau_f  =  length  of  the  virtual  arc 

%  tgo  =  time  to  go  to  intercept 

%  th  =  interceptor  flight  path  angle 

%  thdot  =  rate  of  change  of  interceptor  flight  path  angle 

%  thp  =  first-order  derivative  of  flight  path  angle  wrt  to 

virtual  arc 

%  time  =  optimal  path  time  history  time= [ 0 ; tgo ] 

%  trys  =  collection  of  norms  [X  nX  Xp  nXp  Xpp  nXppp]  (used  for 

plotting) 

%  V  =  velocity 

%  Vdot  =  acceleration 

%  Vp  =  irst-order  derivative  of  velocity  wrt  to  virtual  arc 

%  X, Xp, Xpp, Xppp  =  reference  trajectory  and  its  derivatives  wrt 


tau 

%  xO,xf  =  initial 

%  xdO,xdf  =  initial 
%  xddO,xddf  =  initial 
%  xp, xpp, xppp 
virtual  domain 
%  xt,xdt,xddt 
time=0 


and  final  inteceptor  position 
and  final  interceptor  velocity 
and  final  inteceptor  acceleration 

=  boundary  conditions  (at  0  and  f)  in  the 

=  target  position,  velocity  and  acceleration  at 


global  Re  alphag  path  %  shared  by  SMFlightS,  SMGuidance  & 

SMTra j  ectory 

global  Pos_BR  Pos_SM  Npol  Npt  kpolar  weight90  %  ...  by  SMFlightS  & 

SMTra j  ectory 

global  qq  thdot  psidot  tgo  costs  trys  %  ...  by  SMGuidance  & 

SMTra j  ectory 

%%  Counting  iterations  to  converge 

qq=qq+l; 

tgoold=tgo ; 

%%  Assigning  variables 
tau_f=free (1) ; 
thf  =f ree ( 2 ) ; 
psif  =f ree ( 3 ) ; 

t=BC(10);  %  current  time  in  the  SM  frame  to  compute  it's  stage  (thrust 

and  drug) 

xO=BC  (1 : 3)  ; 
xdO=BC (4 : 6) ; 
xddO=BC  (7 : 9)  ; 

V ( 1 ) =norm (xdO ) ; 
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Vdoti=norm (xddO )  ; 


L(1)=1;%V(1)  ; 
Lpi=0; %Vdoti/L (1) 


xt=BC (11 : 13) ; 
xdt=BC (14 : 16)  ; 
xddt=BC (17 : 19)  ; 


%  Target's  coordinates  at  the  moment  of  detection 
%  Target's  velocities  at  the  moment  of  detection 
%  Target's  accelerations  at  the  moment  of  detection 


V  (Npt )  =2 600-1/3"^  ( tgo-20 )  ;  %  SM  velocity  estimate  at  impact 

Vdotf=-0 . 3 ; %-5 . 9578 ;  %  SM  acceleartion  estimate  at  impact 


L (Npt) =l;%V(Npt) ; 
Lpf=0; %Vdotf/L (Npt) ; 


xf ^xt+xdf^tgo+O  .  5'^xddt'^tgo^2  ; 
xdf= [V (Npt) ^cos (thf ) ^cos (psif ) ; 
V (Npt) ^cos (thf) ^sin (psif) ; 
V  (Npt) ^sin (thf)  ] ; 


thdotf=free (7) ;  psidotf=free (8) ; 

xddf= [Vdotf ^cos (thf) ^cos (psif) -V (Npt) ^cos (thf) ^s in (psif) ^psidotf- . . . 

V  (Npt)  "^sin  (thf)  "^cos  (psif)  "^thdotf; 

Vdotf  "^cos  (thf)  "^sin  (psif)  +V  (Npt)  "^cos  (thf)  "^cos  (psif)  "^psidotf-.  .  . 

V  (Npt)  "^sin  (thf)  "^sin  (psif)  "^thdotf ; 

Vdotf  "^sin  (thf)  +V  (Npt)  "^cos  (thf)  "^thdotf  ]  ; 


%%  Converting  boundary  conditions  ito  the  virtual  domain 
xp0=xd0/L ( 1 ) ; 

xpp0=  (xddO-xdO'^Lpi)  /L  ( 1 )  ^2  ; 
xppp0= [free (4) ; free (5) ; free (6) ] ; 

xpf=xdf /L (Npt) ; 

xppf=  (xddf-xdf "^Lpf )  /L  (Npt)  ^2; 

xpppf= [0; 0; 0]  ; 


%%  Calculating  polynomials'  coefficients  and  reference  trajectories 
dtau  =tau_f/ (Npt-1) ; 
tau  =linspace ( 0 , tau_f , Npt ) ; 
for  i=l:3 
A{i}= [xO (i) ; 

xpO  (i)  ; 
xppO  (i)  ; 
xpppO  (i)  ; 

(-IG^xpppO (i) -4^xpppf (i) ) /tau_f+ (- 
12 0"^ XppO  (i)  TGO'^xppf  (i)  )  / tau_f ^2  +  .  .  . 

(-CGO'^xpf  (i)  -  480"^ XpO  (i)  )  / tau_f ^3+  (840'^xf  (i)  - 
840^x0  (i)  )  /tau_f'"4; 

(  6 0"^ XpppO  (i)  +30'^xpppf  (i)  )  / tau_f  ^2  +  (  GOO'^xppO  (i)  - 
420^xppf  (i)  )  /tau_f'^3+.  .  . 

(2  340^xpf  (i)  +27  00^xp0  (i)  )  /tau_f'^4+  (50  4  0^x0  (i)  - 
5040^xf (i) ) /tau_f^5; 

(-80^xppp0  (i)  -GO^xpppf  (i)  )  /tau_f'^3+  (780^xppf  (i)  - 
900^xpp0  (i)  )  /tau_f ^^4+  .  .  . 
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(-4080'^xpf  (i)  -4320'^xpO  (i)  )  /tau_f^5+  (- 
8400^x0 (i) +8400^xf (i) ) /tau_f^6; 

(35^xppp0 (i) +35^xpppf (i) ) /tau_f^4+ (420^xpp0 (i) - 
420^xppf (i) ) /tau_f^5+. . . 


(2100^xpf (i) +2100^xp0 (i) ) /tau_f "6+ (4200^x0 (i) - 


4200*xf (i) ) /tau_f^7] ; 

Ax{i}  =ciiag  (  [  1,  1, 

1/2, 

1/6, 

1/24, 

1/60, 

1/120, 

1/210 

])*A{i}; 

Axp { i } 

=dlag  (  [  0 ,  1 , 

1, 

1/2, 

1/6, 

1/12, 

1/20, 

1/30 

])*A{i}; 

Axpp { i } 

=dlag  (  [  0 ,  0 , 

1, 

1, 

1/2, 

1/3, 

1/4, 

1/5 

] ) *A{i}; 

Axppp { i } 

=dlag ( [  0 ,  0 ,  0 , 

1, 

1, 

1, 

1, 

1  ])*A{i}; 

Cx (i, : )  =Ax{ i } ( [Npol 
Cxp(i,:)  =Axp { i } ( [Npol : -1 : 2 ] ) ; 

Cxpp (i, : )  =Axpp{ i } ( [Npol : -1 : 3] ) ; 

Cxppp (1, : ) =Axppp{ 1 } ( [Npol : -1 : 4] ) ; 

X  (1,  : )  =polyval (Cx (1,  :) , tau)  ; 

Xp  (1,  : )  =polyval (Cxp (1,  : ) , tau)  ; 

Xpp  (1,  : )  =polyval (Cxpp (1,  : ) , tau)  ; 

Xppp (1, : ) =polyval (Cxppp (1, : ) , tau) ; 

end 

%%  Computing  the  states 
xpl2=Xp(l, :) .-2+Xp(2, :) .^2; 
th  =atan2 (Xp (3, : ) , sqrt (xpl2) ) ; 
psl  =atan2 (Xp (2 ,  : ) , Xp ( 1 ,  : ) )  ; 

thp  =  (Xpp (3,  : )  . ^xpl2- 

Xp(3,  :  )  .^  (Xp(l,  :  )  .^Xpp(l,  :  )+Xp(2,  :  )  .^Xpp(2,  :)))./... 

sqrt (xpl2 ) . / (xpl2+Xp ( 3 , : ) . ^2 ) ; 

pslp  = (Xp (1, : ) . ^Xpp (2, : ) -Xpp (1, : ) . ^Xp (2, : ) ) . /xpl2; 
time ( 1 ) =t ; 

g  =  3.986004418el4/norm(X(l:3,l) ) ^2; 

nx(l)  =  Vdotl/g+sln ( th ( 1 ) ) ; 
ny(l)  =  V  ( 1 ) /g'^cos  ( th  ( 1 )  )  "^psldot  ; 
nz  (1)  =  V(l)/g'^thdot+cos(th(l)); 

[ro, press, temp] =STatmos (norm (X ( : , 1) ) -Re) ; 

[m_l, Sref ] =SMParams3 (time ( 1 ) ) ; 

alphag  (1)  =1 80 /pi sqrt  (ny  (1)  ^2+nz  (1)  ^2)  '^m_l'^g/  (ro'^V  (1)  ^2'^Sref/2)  /13; 
for  j  =2 : Npt ; 

g=3.98  6004  418el4/norm(X(l:3,  j)  )  ^2; 

If  norm (X (:, j )) -Re<8 6000 ,  [ ro, press , temp] =STatmos (norm (X (:, j )) - 

Re)  ; 

else  [ro, press, temp] =STatmos (86000) ; 

end 

MV=V( j-1) /sqrt (1.402^287. 053^ temp) ; 

CDtable=SMDrag (MV) ; 

If  tlme(j-l)<20  Thrust=95300 ;  CD0=CDtable ( 1 ) ; 

else  Thrust=0;  CD0=CDtable (2 ) ; 

end 

[m_l, Sref] =SMParams3 (time ( j -1 ) ) ; 
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CL=  ( 2  "^sqrt  (ny  (j-1)  ^2+nz  (j-1)  ^2)  /  ( ro'^V  ( j  -1 )  ^2  "^Sref )  ; 

CD  =  CDO  +  kpolar^CL^2; 

Drag^rc^V  (j-1)  ^2'^CD'^Sref/2; 
nx ( j ) = (Thrust-Drag) /m_i/g; 

V(  j  )  =V(  j-1)  +g^  (nx(j-l)-sin(th(j-l)  )  )/L(j-l)  ^dtau; 

ddist^sqrt ( (X(l,j)-X(l,j-l))-2+(X(2,j)-X(2,j-l))-2+(X(3,j)-X(3,j- 
1))"2); 

dtime=2^ddist/ (V ( j ) +V ( j -1 ) ) ; 

L ( j ) =dtau/ dtime ; 

ny ( j ) =V ( j ) /g^cos (th ( j ) ) ^psip ( j ) ( j ) ; 
nz ( j ) =V ( j ) /g^thp ( j ) ( j ) +cos (th ( j ) ) ; 

alphag  ( j  )  ^ISO/pi'^sqrt  (ny  (j)  ^2+nz  (j)  ^2)  '^m_i'^g/  (rc^V  ( j  )  ^2'^Sref/2)  /13; 
time ( j ) =time ( j -1 ) +dtime ; 

end 

tgo=time (end) -time (1) ; 

for  i=l : Npt 

nX ( i ) =norm (X ( 1 : 3 ,  i )  )  ; 
nXp ( i ) =norm (Xp ( 1 : 3 , i ) ) ; 
nXpp (i) =norm (Xpp (1 : 3, i) ) ; 

end 

trys=[X'  nX '  Xp '  nXp '  Xpp'  nXpp ' ] ; 

path= [time '  X(l:3,:)'  V  th'  psi'  nx'  ny '  nz '  tau '  L ' ] ; 

%%  Computing  cost  and  penalties 
V_f=V(end) ; 

xdt=BC (14 : 16) +BC (17 : 19) ^tgo; 
xdf=  [V_f "^cos  (thf )  "^cos  (psif )  ; 

V_f "^cos  (thf)  "^sin  (psif)  ; 

V_f^sin (thf)  ]  ; 

nmax=40+ (40-10) / (0-50000) ^ (X (3, : ) -Re) ; 
if  nmax<l,  nmax=l;  end 

J= [ tgo; 

abs (dot (xdf , xdt) ) /norm (xdf ) /norm (xdt) ] ; 

Py=max ( [ 0 , abs (ny) -nmax] ) ; 

Pz=max ( [ 0 , abs (nz ) -nmax] ) ; 

Dtgo=tgoold-tgo ; 

Altf  ine=max  (  [  0 ,  X  ( 3 ,  end)  -Re-60000  ]  )  .  '^P+min  (  [  0 ,  X  ( 3 ,  end)  -Re-9000  ]  )  .^2; 
cost=norm  {[1 ,  weight  90  ]  J)  +10"^  (Py^2  +  Pz^2 )  TlOO'^Dtgo^P  +  O  .  I'^Altf  ine; 
costs (qq, : ) = [cost; tau_f ; J; Py; Pz ; Dtgo ; Altf ine ] ; 

%%  Animating  iterations 


f igure ( 100 ) , %  set (gcf ,  'Color '  ,  ' w ' ) ; 
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subplot (3, 6,  [1  2  7  8] ) 

plots (Pos_BR(l:300, 1) /10"3 , Pos_BR  ( 1 : 300 , 2 ) /10"3,  (Pos_BR (1 : 300, 3) - 
Re) 710^3,  .  .  . 

' - . r ' , ' LineWidth ' , 2 ) 

hold  on,  ylim([0  60]);  %axis  equal 

plots (Pos_BR( 14 0,1) 710^3, Pos_BR( 14 0,2) 710^3, (Pos_BR (140, 3) -Re) 710^3, . . . 

''^g'  ,  'LineWidth'  ,2) 

plots (Pos_SM(l : 19, 1) 710^3, Pos_SM (1 : 19, 2) 710^3, ( Pos_SM ( 1 : 1 9 , 3 ) - 
Re) 710^3,  .  .  . 

' . k ' , ' LineWidth ' , 2 ) 

plots (X(1,:)710^3,X(2,:)710^3, (X(3,:) -Re) 710^3, 'Linewidth' , 3) ;  grid  on 
plots (xf (1) 710^3, xf (2) 710^3, (xf (3) -Re) 710^3, 'pk' , 'MarkerSize' ,11) 
plots (X(l,end) 710^3, X (2 , end) 710^3,  (X(3,end) - 
Re) 710^3,  'pr',  ' MarkerSize ' , 9)  ; 
view (-102 , 8 ) 

hl=legend ( ' BM  tra j ectory ' , ' BM  detection ',' Unguided  ascend',... 

'Guided  flight ',' Predicted  intercept  point ',' Actual  intercept 
point ' , . . . 

' Location ' , ' North ' ) ;  set (hi , ' Font Size ' , 7 ) 
hold  off 
%{ 

axis([0  l.le5  -le4  1 . Ie5  0  7e4]7l0"3) 

xlabel (' Northing,  x  ( km)  '), ylabel (' Westing,  y  ( km)  '), zlabel (' Altitude 
(km)  ') 

subplot (3, 6, [13  14] ) 

plot (path ( : , 1 ) , path ( : , 9 ) , ' --b ' , ' Linewidth ' , 2 ) 
hold  on;  grid  on 

plot  (path ( : , 1 ) , path (:,10),  '-.g',  ' Linewidth ' , 2 ) 
plot (path ( : , 1 ) , nmax ( : ) , ' r ' , ' Linewidth ' , 2 ) 
plot (path ( : , 1 ) , -nmax ( : ) , ' r ' , ' Linewidth ' , 2 ) 
axis ([10  time (Npt)  -45  45]); 

xlabel ('Time  (s)'),  ylabel ('Load  Factor  (g)  '  ) 
leg2=legend ( ' n_y ' , ' n_z ' , ' Dynamic 
constraints ' , ' Location ' , ' Best ' ) ; 
set(leg2, ' FontSize ' , 7 ) 
hold  off 
subplot (3,6,3) 

plot (time, X ( 1 ,:) 710^3 ,' Linewidth ', 2 ) ;  grid  on 
axis([10  time(Npt)  Ie47l0"3  I.le57l0"3]) 
title ( ' x_l  ( km) ' ) 
subplot ( 3 , 6 ,  9 ) 

plot (time, X (2 ,:) 710^3,  ' Linewidth ', 2 ) ;  grid  on 
axis([10  time(Npt)  -Ie47l0'"3  1 .  Ie57l0'"3  ]  ) 
title ( ' x_2  ( km)  '  ) 

subplot ( 3 , 6 , 15 ) 

plot (time, (X ( 3 , : ) -Re) 710^3 , ' Linewidth ' , 2 ) ;  grid  on 
axis  ([10  time  (Npt)  0  7e47l0"^3]) 
title ('x_3  (km)'),  xlabel ('Time  (s)') 

subplot (3,6,4) 

plot  (time,  Xp  ( 1 ,:)  710^^3 ,' Linewidth ',  2 )  ;  grid  on 
axis ([10  time (Npt)  -45  45]);  axis  'auto  y' 
title('x_l'  '  710'"3') 
subplot (3,6,10) 
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plot  (time,  Xp  (2 ,  :  ) /10"^3,  '  Linewidth  '  ,  2 )  ;  grid  on 
axis ([10  time (Npt)  -45  45]);  axis  'auto  y' 
title('x_2' '  /10"3') 
subplot (3,6,16) 

plot  ( time,  Xp  ( 3 ,  :  ) /10'^3 ,  '  Linewidth  '  ,  2 )  ;  grid  on 
axis ([10  time (Npt)  -45  45]);  axis  'auto  y' 
title  ('x_3''  710^3'),  xlabeK'Time  (s)') 
subplot (3,6,5) 

plot  ( time,  Xpp  ( 1 ,  :  ) /10'^2 ,  '  Linewidth  '  ,  2 )  ;  grid  on 
axis ([10  time (Npt)  -45  45]);  axis  'auto  y' 
title ( 'x_l ' ' ' '  710^2 ' ) 

subplot (3,6,11) 

plot  (time,  Xpp  (2 ,:)  710^^2 ,' Linewidth ',  2 )  ;  grid  on 
axis ([10  time (Npt)  -45  45]);  axis  'auto  y' 
title ( 'x_2 '  '  '  '  710^2  '  ) 

subplot  ( 3 , 6 , 17 ) 

plot  ( time,  Xpp  ( 3 ,  :  )  710"^2 ,  '  Linewidth  '  ,  2 )  ;  grid  on 
axis ([10  time (Npt)  -45  45]);  axis  'auto  y' 
title  ( 'x_3'  ''  '  710^2'),  xlabelCTime  (s)') 
subplot  (3, 6, 6) 
hold  on 

plot (qq, time (Npt ) -time ( 1 ),' +c ',' Linewidth ', 1 ) ;  grid  on 
hold  off 

axis([l  200  40  60]);  axis  'auto  y' 
title (' Time-to-go  (s)') 
subplot (3, 6, [12  18] ) 
hold  on 

plot (qq, cost ,' +m ',' Linewidth ',  1 )  ;  grid  on 
axis([l  200  0  100]);  axis  'auto  y' 
hold  off 

title (' Performance  index'),  xlabel (' iteration ' ) 

%} 

return 


D.  COMMON  FUNCTION  -  STANDARD  ATMOSPHERE 

1.  STatmos.m 


function  [Density,  Pressure,  Temperature] =STatmos (alt ) 

%  Calculation  of  the  1976  standard  atmosphere  up  to  86  km 
%  Code  source:  http:77www.pdas.com7atmos.htm 

%  Run  ezplot ( ' STatmos ' , [ 0 , 8 6000 ] )  to  see  the  plot  of  density  vs 
altitude 

o 

o 

%  Author:  Yakimenko,  Oleg  A. 

%  Date:  September,  27  2005 

%  E-mail:  oayakime@nps.edu 

o 

o 

alt=alt7l000 ;  %  Convert  altitude  from  m  to  km 

%%  -  initialize  values  for  1976  atmosphere 
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REARTH=6369 . 0 ;  %  Earth  radius  (km),  depends  on  Latitude 

GMR=34 . 163195;  %  Gas  Constant 

htab=[0.0,  11.0,  20.0,  32.0,  47.0,  51.0,  71.0,  84.852];  %  Geometric  alt 
ttab= [288 . 15,  216.65,  216.65,  228.65,  270.65,...  %  Temperature 

270.65,  214.65,  186.946]; 

ptab=[1.0,  2.233611E-1,  5.403295E-2,  8 . 5666784E-3, . . .  %  Relative  pres 

1 . 0  94  5  601E-3,  6 . 60  63531E-4,  3 . 904  6834E-5,  3 . 6850 IE- 6] ; 
gtab=[-6.5,  0.0,  1.0,  2.8,  0.0,  -2.8,  -2.0,  0.0];  %  Temp  gradient 

P0=101325.0;  Ro0=1.225; 

%% -  Convert  geometric  to  geopotential  altitude 

if  alt>250 

alt  =  100 

end 

h=alt^REARTH/ (alt+REARTH) ; 

%%  -  Binary  search  for  altitude  interval 

i  =  1; 

j  =  8; 

while  j  >  i+1 

k=fix ( (i+j ) /2) ; 
if  h<htab ( k) ; 

j  =  k; 
else 
i  =  k; 
end 

end 

%%  -  Calculate  local  temperature 

tgrad  =  gtab  (i) ; 

tbase  =  ttab  (i) ; 

deltah  =  h  -  htab(i); 

tlocal  =  tbase  +  tgrad^deltah; 

theta  =  tlocal/ttab ( 1 ) ; 

%%  -  Calculate  local  pressure 

if  (tgrad  ==  0.0) 

delta  =  ptab  (i) ^exp (-GMR^deltah/tbase)  ;  %  Isothermal  layers 

else 

delta  =  ptab  (i)  ^  (tbase/tlocal) (GMR/tgrad)  ;  %  Non-isothermal 

layers 
end 

%%  -  Calculate  local  density 

sigma  =  delta/theta; 

%%  -  Current  atmosphere  parameters  corresponding  to  Altitude  alt 

Temperature=tlocal ;  Density^RoO'^sigma;  Pressure^PO"^  delta; 
return 
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APPENDIX  E:  DETAILED  DESCRIPTION  OF  TRAJECTORY 
SHAPING  GUIDANCE  (REPRODUCED  FROM  [10]) 


A  method  that  overcomes  fatal  characteristics  of  modem  guidance  laws  is  the  key 
driving  factor  in  the  development  of  this  advanced  guidance  law.  The  guidance  law  will 
determine  the  near-optimal  flight  path  from  the  interceptor  position  to  a  predicted  target 
position  for  the  interceptor  to  follow  to  intercept,  and  then  derive  the  set  of  control 
commands  necessary  to  execute  that  flight.  This  section  will  describe  a  method  for 
deriving  that  trajectory  using  calculus  of  variations  based  on  three  cues:  high-order 
polynomials  as  a  reference  function  for  the  flight  path,  a  preset  thmst  history  as  one  of 
the  controls,  and  a  few  optimization  parameters.  The  trajectory  optimization  problem  is 
then  converted  into  a  nonlinear  programming  problem  and  solved  numerically. 

A.  PROBLEM  STATEMENT 


Among  all  admissible  trajectories. 


that  satisfy: 

zit)  =  {z^it),z^{t),...,zXt)Y 

S  =  {z{t)GZ’- 

(4.A.1) 

1. 

The  system  of  differential  equations  (dynamic  constraints): 

Zi  =  fiit,z,u,a),  i  =  l,r 

where  the  vector  of  controls  is 

(4.A.2) 

u(t)  =  ,  m<r,  ugU"'<^E"’  and  the 

missile  parameters  is  a  =  {a^,a2,...,ap),  a  e  E^ ; 

vector  of 

2. 

The  initial  conditions: 

z(to)eSo,  S,={z,eZ^^EY 

(4.A.3) 

uit,)eR„  R,={u,eU"‘c^E”} 

and  the  final  conditions: 

(4.A.4) 

z(t,)eS,.  S,^{z,^Z‘  G,{z{t,))  =  0,j  =  ri} 

(4.A.5) 

m 

II 

m 

s 

n 

(4.A.6) 

3. 

The  constraints  imposed  on  the  state  space 
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fj{t,  z)  =  (t,  z),  ^2 (^,  z), ...,  (t,  z)}^  >  0 

(4.A.7) 

on  the  controls 

i(t,Z,u)  =  {^^(t,Z,u),^2(^,Z,u),...,<^^(t,Z,u)f  >6 

(4.A.8) 

and  on  the  controls  derivatives 

(4.A.9) 

Find  the  optimal  trajectory,  z^^^{t) ,  that  minimizes  the  integral  function 

‘f 

J  =  K{x^,Xj)  + ^  Lit,z,u)dt  (4.A.10) 

^0 

and  the  corresponding  optimal  controls,  (t) ,  where  K,  L  are  defined  functions. 

B.  CALCULUS  OF  VARIATIONS 

Calculus  of  variations  deals  with  functions  of  functions,  termed  functional, 
instead  of  functions  of  some  variable  or  variables  as  in  ordinary  calculus.  Specific 
interest  is  in  the  externals  of  these  functionals  -  those  making  the  functional  attain  a 
maximum  or  minimum  value  [Ref  15].  There  are  two  broad  categorizations  of  methods 
to  solve  these  problems,  indirect  methods  and  direct  methods. 

Indirect  methods  resolve  the  problem  into  a  differential  equation,  usually  via  the 
difference  of  a  series  of  equations  of  motion  and  thus  solve  the  general  theory  of  partial 
differential  equations.  This  method  does  not  assume  anything  about  the  solution,  leading 
to  a  very  complete  and  perfectly  precise  answer;  however,  one  must  integrate  the 
resultant  equation  to  derive  the  extremals.  The  precision  greatly  increases  the 
computational  complexity,  and  therefore  the  time  required  to  solve,  for  negligible  gain  in 
optimality.  Further,  these  differential  equations  are  difficult  to  integrate  except  in  the 
simplest  of  cases,  and  nearly  impossible  to  program  [Ref  19].  This  approach  is  further 
complicated  by  the  need  to  solve  the  problem  in  a  specified  fixed  region  instead  of  in  the 
small  neighborhood  of  some  point.  These  difficulties  can  be  overcome  by  using  direct 
methods,  which  do  not  reduce  the  variational  problems  to  ones  involving  differential 
equations. 
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The  fundamental  idea  of  direct  methods  is  to  consider  a  variational  problem  as  a 
limit  problem  of  the  extreme  of  a  function  of  a  finite  number  of  variables  that  can  be 
solved  by  numerical  methods.  Two  basic  direct  methods  are  the  Rayleigh-Ritz  and 
Galerkin  methods,  which  assume  the  solution  to  be  an  unknown  function  but  of  a  certain 
form  containing  a  set  of  unknown  coefficients  (themselves  functions  of  the  boundary 
conditions),  which  are  then  found  by  minimization.  The  practical  result  is  that  the 
problem  has  been  reduced  from  calculus  to  algebra,  but  at  the  cost  of  a  significant 
increase  in  the  number  of  simultaneous  equations  to  be  solved.  It  is  for  this  reason  that 
little  work  was  done  in  this  field  until  the  advent  of  computer  technology.  The  possible 
solution  set  is  restricted  to  a  smaller  space  than  the  original  equation  because  of  the  initial 
assumption  regarding  the  solution,  but  the  resultant  problem  can  be  programmed  and 
quickly  solved  using  computers;  however,  the  solutions  are  only  an  approximation  of  the 
original  solution.  Therefore  these  solutions  can  only  be  regarded  as  near-optimal 
solutions. 


Professor  Taranenko  [Ref  19]  first  applied  the  ideas  of  direct  methods  and  the 
combination  of  Ritz  and  Galerkin  methods  to  the  problems  of  flight  dynamics  by 
identifying  a  reference  function  for  the  flight  vehicle’s  motion  and  velocity 


X.  =x,.o+(x.-x,.o) 


r-n 


Tf  Tq 


+  <I),(r),  1  =  1,4 


(2.1) 


where  Xj,X2,X3,are  the  Cartesian  coordinates  for  the  flight  path,  x^is  the  velocity,  and 
0.(r)  is  a  continuous,  unequivocal,  differentiable  function  satisfying  the  boundary 


conditions  <l),.(ro)  =  O.(ry^)  =  0  .  Any  function  that  satisfies  those  conditions  would  be 

acceptable  for  use.  Taranenko  called  r  a  virtual  arc.  It  is  this  critical  variable  that  allows 
the  separation  of  the  spatial  trajectory  from  the  velocity  and  thus  optimize  one  or  the 
other,  or  both  independently.  The  specific  task  determines  the  appropriate  choice  of  t  , 
but  in  general  any  continuous,  monatomic  function  is  acceptable:  time,  path,  energy,  etc. 
The  remaining  state  parameters  and  flight  controls  are  then  determined  by  solving  the 
inverse  problem  of  flight  dynamics.  Rather  than  starting  with  the  control  time  histories 
and  integrating  them  to  determine  the  flight  path,  this  method  starts  with  a  flight  path  and 
determines  the  control  time  histories  necessary  to  create  it. 
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In  order  to  implement  a  solution  to  this  problem  that  can  be  solved  in  real  time,  a 
further  restriction  must  be  made.  The  continuous  problem  must  be  discretized  to  reduce 
the  infinite  variational  problem  to  one  of  optimization  of  few  parameters  at  numerous 
sampling  points.  This  allows  for  the  optimization  of  the  planar  trajectory  of  the  missile  at 
several  points  along  the  path  by  presetting  the  state  variables  and  one  of  the  controls’ 
time  histories,  and  then  solving  the  inverse  flight  dynamics  problem. 

These  methods  have  all  previously  been  used  for  off-line  optimization  of 
trajectories,  but  none  has  yet  been  applied  to  the  real-time  onboard  optimization  of  a 
missile  flight  path.  Yakimenko  detailed  a  method  called  a  Direct  Method  for  Rapid 
Prototyping,  which  he  applied  to  short  term  spatial  trajectories  of  aircraft  maneuvers 
using  fixed  boundary  points  [Ref  19].  The  developed  program  presented  here  uses  similar 
numerical  method  to  provide  a  near-optimal  spatial  trajectory  that  is  completely  defined 
by  a  few  optimization  parameters,  but  as  will  be  discussed  shortly,  has  fluid  final 
boundary  conditions. 

Though  the  method  artificially  limits  the  possible  trajectory  variations,  it  does 
guarantee  the  following  [Ref  19]: 

1 .  The  boundary  conditions  are  satisfied  a  priori, 

2.  The  control  commands  are  physically  realizable  and  smooth, 

3.  Only  a  few  variable  parameters  are  used,  thus  ensuring  that  the  iterative 
process  converges  well, 

4.  The  near-optimal  solution  is  very  close  to  the  optimal  one. 


A  simple  two  dimensional  variation  program  of  a  similar  7*  order  system  will 
demonstrate  how  the  direct  method  varies  the  flight  path  according  to  the  boundary 
conditions.  The  boundary  conditions  are 


o 

II 

o 

o 

II 

o 

=  1 

X2/  - 1 

x[q  =  0.2 

^20  1 

Xj'/  = 

0.1 

X2y  =  -1 

II 

p 

II 

P 

ff 

X,f  = 

0.1 

x"^  =0.1 

x"\q  =  var 

X"'20  =0.1 

x'"  = 

0.1 

x"'2^  =0.1 
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where  the  third  derivative  of  the  initial  Vj  condition  has  been  set  to  vary  according  to 

=  {-0.4; -0.1;  0.2;  0.5} 

and  the  length  of  the  virtual  arc  varies  according  to 

=1,2,. ..,10 

The  resulting  set  of  paths  is  shown  in  figure  26  and  the  first  derivative  of  the  path  is 
shown  in  figure  27.  It  is  important  to  note  that  the  first  derivative  is  not  “velocity”  but  is 
instead  the  “rate  of  change  of  the  path”.  It  is  proportional  but  not  equal  to  velocity, 
specifically  because  of  the  virtual  variable  r  as  discussed  previously.  Each  line 
represents  a  different  choice  of  ,  showing  that  by  varying  that  value  the  length  of  the 

path  can  change  drastically.  The  algorithm  developed  here  will  use  this  technique  to 
derive  the  optimal  flight  path. 


dVdT^|g=-0.4 


T^=var 

0  0.5  1  1.5  2 

d\ldx\=0.2 


d\/dT^|Q=-0.1 


d\/dT^|g=0.5 


Figure  37.  Variation  of  Path  with  Tj-  and  x}"  (after  [Ref  18]) 
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d\/dT^|Q=-0.1 


3 

2 

1 

0 
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d3x/dT3|Q=0.5 


T 

Figure  38.  Variation  of  First  Derivative  of  Path  with  and  Xj" 


C.  PROGRAM  DEVELOPMENT 

1.  Boundary  Conditions 

Using  equations  (3.2)  and  (3.3),  in  addition  to  the  controls  n^,ny,n^,  it  is  possible 
to  construct  the  vector  of  state  variables  z  =  {xj,X2,X3,F,^,'P}^  and  the  vector  of  controls 
u  =  {n^,ny,n^Y .  The  model  for  thrust,  drag,  and  missile  characteristics  is  the  same  as 
previously  used  in  Chapter  3. 

The  beginning  assumption  is  that  the  following  data  is  known  by  the  interceptor’s 
onboard  computer: 
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Interceptor 

Target 

Body  Frame 

Xi(^o)  =  ^10 
^2  (^0  )  “  ^20 
^3  (^0  )  “  -^30 

^(^o)  “  ^0  -^10 

^(^o)  “  ^0  ^  ^20  ^10 

^(^o)  ~  ^0  -^30  ^-^20 

^y(to)  =  ^yO  ^30 

Xi(t/)  =  Xi^ 

■^2  (^/  )  ~  -^2/ 

■^3  (^/  )  ~  -^3/ 

V{t^)  =  Vf  'I  x,^ 

'P(t^)  =  'P,  k/ 

)  ~  _ .  -^3/ 

Earth  Centered 
Inertial 

xf^(t),  xf^it),  xl^(t) 

X™(t),  xfCt),  X3""(t) 

x“(t),  xf  (t),  xl^(t) 

xf  (t),  xf  (t),  xf  (t) 
xf"(t),  xf  (t),  X3^"(t) 

Table  9.  Interceptor  Known  Data 


The  target  data  will  be  used  to  determine  the  final  boundary  conditions  of  the 
interceptor  missile  according  to  the  time  until  intercept,  At : 

Position: 


SM 


BR 


x^y  =  Xjo '  +  Xjo  At 


SM 


BR 


^^'2.  •X^2o  I  *^20 


.,SM 


BR 


j"  *^20  ^  *^20 


Heading  and  Flight  Path  Angle: 


gSM  ^  :i  QBR  ^  ^  diVCig 


n 


BR2  ,  • BR2 

\  *^2 


•BR 


wT 


Xr, 


DU  ^  DD  ^ 

y/j-  ,  where  y/j-  =arctg- 


X, 


(4.C.1) 


(4.C.2) 


which,  when  combined  with  reasonable  estimates  of  the  final  values  of  the  velocity,  V, 
and  the  time  rate  of  change  of  velocity,  V ,  heading  ( 'P  =  0 ),  and  flight  path  angle 
(0  =  0),  yields  the  final  conditions: 
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(4.C.3) 


xf^  =  Vj-  cos  cos  \j/ f 

■^2/  =  ¥f 

xf!  =V^  sine f 

and 

x™  =  cos  Of  cos  -  OfVf  sin  Oj  cos  ^  ^  ^  fVf  cos  Of  sin 

x™  =  Vf  cos  Of  sin  f  -  OfVf  sin  Of  sin  +  T^F^  cos  Of  cos  f  (4.C.4) 
x™  =  Vf  sin  Of  +  OfVf  cos  Of 

In  order  to  ensure  a  smooth  flight  path  at  the  initial  and  final  points,  an  additional 
constraint  of 


■■■SM  _  ■■■SM 
■^U  ~  -^1/ 


X. 


.■SM 


=  X, 


X, 


■SM 


=  X 


.■SM 

■2/ 

SM 

•3/ 


=  0 
=  0 
=  0 


will  be  imposed  on  the  system  at  the  initial  and  final  eonditions. 


(4.C.5) 


2.  Separating  and  Recombining  Space  and  Time 


As  discussed  previously,  in  order  to  independently  optimize  the  spatial  trajeetory 
and  the  veloeity,  the  reference  funetion  will  be  derived  as  a  function  of  r  .  The  boundary 
eonditions  eannot,  therefore,  be  defined  as  funetions  of  time  derivatives  as  in  equations 
(4.C.3)  -  (4.C.5).  A  eonnection  between  the  spatial  and  time  domains  must  therefore  be 
introduced,  T ,  which  is  defined  as 


(4.C.6) 


and  is  termed  the  virtual  speed  [Ref  19],  This  allows  for  the  independent  variation  of  the 
speed  profile  along  the  same  paths  aceording  to  any  other  convenient  referenee.  In  this 
ease  the  known  thrust  profile,  will  be  utilized  by  integrating  the  third  equation  of 
(2.B.3)  and  applying  the  virtual  speed 


jM,  ^  m  dr  g{n^- sin  0) 

V  {T)  =  g{n^-smO)—  = - - 

at  a{t) 


(4.C.7) 
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Further  use  of  the  virtual  speed  allows  the  recalculation  of  the  initial  and  final  boundary 
conditions,  transforming  them  from  the  time  frame  to  the  spatial  frame.  For  this,  the 
obvious  relations 


Xi{T)  =  ^^  =  x'.{T)X{T) 

cIt  at 


x,.(r)  = 


d{x'.{T)X{T))  dr  „.2 


(4.C.8) 


=  x”a  +  i  =  1, 2, 3 


dr  dt 

which  when  rearranged  defines  the  first  and  second  derivatives  of  the  missile  coordinates 
as 

x;  =  r‘x,.  x;'=r'[x,.-x.A']  / =1,2,3  (4.C.9) 

Using  the  values  of  A,  and  X'  defined  as 


X  =  V  X'=VV  X  =V  X'  =  V  V 

/u  'O’  '‘"0  'O'O  '‘'f  '  f'  t 


-1 


(4.C.10) 


3.  Reference  Trajectory 


The  knowledge  of  the  initial  and  final  position  plus  the  initial  and  final  conditions 
of  the  first  and  second  time  derivates  allows  for  the  construction  of  a  7*-order 
polynomial  (the  maximum  orders  of  the  time  derivatives  of  the  missile  coordinates  at  the 
initial  and  final  points  plus  one)  to  describe  the  reference  function  of  the  aircraft 
coordinates  x.  (/  =  1, 2, 3) .  The  following  are  introduced  as  the  reference  functions  [Ref 


19]: 


Or  written  another  way. 


(max(l,A:-2))!r^ 


ik 


k\ 


k=0 


k  =  l 
5 


a.,T 


k-2 


■(r> = s 

k=2 

■;Tr)  =  2(t-2)a„r 


1-3 


1=3 


(4.C.11) 
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—  Cl-^  +  Cl^^T  + 

x"i{T)  =  ai^+ai^T  +  -anT^  +-ai^T^ 

x'.(r)  =  a,.;  +a;2r +  -a;3r^  +7«;4^^  +7:7^/5^''  +T7:«,7^^ 

2  6  12  20  30 

,  ^  1  2  1  3  1  4  1  5  1 

x,(r)  =  a,o+fl,ir  +  -fl,2r  +-a,3r 


(4.C.12) 


fl.,r  + 


210 


a.^T 


The  coefficients  can  be  determined  by  solving  the  equations  simultaneously, 


fl 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

2  V 

6  V 

J_7-4 

24  ‘■f 

60  .f 

120  V 

2T0 

0 

1 

2  V 

6  V 

J_7-4 

12  V 

20  V 

30  V 

0 

0 

1 

12-2 

2  V 

if 

i^/ 

5 

lo 

0 

0 

1 

f 

a-o 

a., 


a.2 


a, 4 


V6 


^  X  ^ 
■^iO 


^,0 

ff 

^,0 


X 


X,, 


X, 

x; 

x" 


(4.C.13) 


V-"  ifj 

Finally,  by  substituting  the  corresponding  values  of  x^q,x\^,x”^  (z  =  1,2,3)  for 
To  =  0,  and  x.j-,x'.j-,x'l^  (i  =  1,2,3) for  t  =  Tj- (where  Tj- is  the  first  optimization  parameter, 

the  virtual  arc),  results  in  a  set  of  24  linear  algebraic  equations  for  21  unknown 
coefficients  a.,^  (z  =  1,2,3,  A:  =  0,1,.., 7) 


^iO  ^iO 


a.j  x.Q 


a.3  =  X 


-4x!'  +  16x'q  60x,!^ -120x'q  -360x7 -480x'f,  840x..  -  840x^0 

«,-4  = - ^ - + - S - + - S - 


V 


V 


V 


SOx^  +  bOx';;  -420x" +600x7  2340x7  -  2700x'  -5040x,,  +  5040x,„ 

a,.,  = - ‘:L_ - ^  + - £_ - ^  + - 'L_ - ^  + - 'L_ - ^(4.C.14) 


^i6  = 


V 


V 


V 


-60x^-80x;'  780x^-9004  -4080x;^ -4320x;o  8400x^^  -  8400x.o 

7  A  /Z  if 


V 


V 


V 


SSx^'  +  354'  -420x;:  +  420x;'o  2100x;o  +  2100xL  -4200x.,  +  4200X.O 

S —  + - S - + - ^ ^  + - S - 


V 


V 


V 


'f 


no 


4.  Inverse  Dynamics 

The  numerical  solution  develops  a  reference  trajectory  over  a  fixed  set  of  points 
equidistantly  placed  along  the  virtual  arc.  The  virtual  interval  is 


Ar  =  — — 
A^-1 


(4.C.15) 


and  the  corresponding  time  interval  is 


f  3 


=  2 


V  1  _ y 


(4.C.16) 


where  Fj  and  is  determined  by  integrating  equation  (4.C.7).  N  is  any  convenient 
number,  chosen  to  be  100  in  this  program. 

The  value  of  Fj  also  allows  for  the  determination  of  both  6*  and  'P  by  rearranging 
equations  (2.B.2)  and  substituting  the  virtual  velocity 

f 

X.  . 

6.  =  arctg- 

Jxl.f+xf.." 

(4.C.17) 


X. 


=arctg— 


i\j 


The  controls  and  are  found  by  rearranging  equations  (2.B.3)  to  be 

F.  . 

n^.j  =—^j  CO?,  Oj 

S 

n  =^{6  +  cos  6) 

g 

where  the  angular  derivatives  are  determined  as 

0  =  COG^  9  ~  -^3 (-^1'+  ^2) 

'  (x;^+x;^f^ 

-  \T/  ^1^2  ~  ^1^2  ) 


'P,  =cos"T'- 


X, 


,.'2 

1 


(4.C.18) 


(4.C.19) 
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5. 


Cost  and  Penalty  Functions 


Finally,  the  calculation  of  the  flight  path  results  in  a  set  of  functions  that  must  be 
minimized,  which  occurs  through  a  Cost  Function  (CF)  and  a  Penalty  Function  (PF). 
These  functions  must  be  carefully  chosen  to  ensure  that  the  optimal  path  is  truly  feasible, 
desirable,  and  obtainable.  A  simple  example  of  a  CF  is  J  =  1^ ,  the  minimal  time 

problem,  ox  J  =  |M(t)| ,  the  minimum  fuel  problem,  though  there  is  no  limit  to  the  number 
or  variation  of  the  Cost  Function. 


In  this  case,  the  CF  was  chosen  to  optimize  three  properties  simultaneously; 
minimize  the  length  of  the  virtual  arc,  minimize  the  time  to  intercept,  and 

maximize  the  impact  angle  of  the  interception.  The  CF  for  this  program  is  written  as 


J  =  W^k^T j  +  W2k2tgg  +  W^kj 


-V 


SM 


(4.C.20) 


Each  item  must  be  scaled  appropriately  using  the  scaling  factors  k^,k2,k 2,  so  that 


they  are  all  roughly  equivalent  when  optimized,  e.g.  the  anticipated  intercept  time  is 
counted  in  tens  of  seconds  while  the  cosine  of  the  impact  angle  will  vary  from  zero  to 
one.  Failure  to  weight  them  properly  will  skew  the  results  of  the  cost  function. 
Additionally,  through  the  weighting  functions  Wj,W2,W3,a  trade-off  analysis  can  be 
conducted  and  variables  can  be  included  or  excluded  as  desired. 


The  first  two  variables,  t  and  are  necessary  to  ensure  the  system’s  optimal 

solution  is  actually  physically  realizable.  Without  them,  the  system  will  continue  to 
optimize  the  intercept  well  beyond  the  capabilities  of  the  missile  or  even  physical  reality. 
One  example  is  that,  in  a  purely  mathematical  sense,  there  is  no  problem  with  a  negative 
velocity  (yet  in  the  physical  world  that  makes  no  sense)  and  the  program  might  try  a 
program  that  would  intercept  after  the  missile  velocity  has  gone  past  zero  and  into 
negative  numbers  (of  course,  the  missile  would  stop  flying  long  before  it  even  reaches 
zero).  One  might  be  tempted  to  include  a  myriad  of  parameters  to  cover  all  possible 
eventualities;  however,  including  these  two  parameters  sufficiently  accounts  for  nearly 
every  physical  limitation  and  eliminates  the  need  for  having  a  long  list  of  cost  variables. 
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The  third  variable  in  the  CF  was  chosen  in  order  to  maximize  the  angle  of  impact 
upon  interception.  This  reflects  the  need,  described  in  Chapter  1,  to  maximize  the  kinetic 
energy  in  order  to  ensure  the  interceptor  disables  the  target.  The  cost  is  calculated  using 
a  simple  dot  product  relationship 


cos^  = 


r  >  /  SM  BR  \ 

dot{Xf  ,  Xj-  ) 


ZSM 

Zbr 

Xf 

(4.C.21) 


which  will  have  a  minimum  value  of  zero  when  the  impact  angle  is  maximum. 


The  PF  is  chosen  to  ensure  that  certain  conditions  are  not  violated  or  exceeded, 
such  as  physical  limitations.  In  this  case,  the  PF  is  on  the  maximum  acceleration  in  the  y 
and  z  direction.  Zarchan  showed  that  acceleration  capability  is  dependent  on  altitude  and 
speed  [Ref  21].  The  PF  has  been  set  up  to  reflect  this  by  varying  between  40  g’s  at  sea 
level  to  10  g’s  at  50,000  ft. 


These  are  not  the  only  CF  or  PF  variables  that  can  be  included.  The  choice  of 
variables  to  include  is  situation  dependant  and  may  be  modified  to  meet  whatever  the 
needs  of  the  situation  demand. 
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